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ABS'TRACT 


New concepts in aerospace travel have renewed interest 
in modeling and compensating for the effects of upper 
atmospheric drag. In particular, the SDI constellation 
requires strict orbital element maintenance. This thesis is 
a qualitative verification of a propellant longevity model 
for low-altitude earth orbit satellites doing intrack micro- 
thrusting to overcome atmospheric drag. The original model 
was developed at the Naval Surface Weapons Center. 
Pertinent orbital mechanics and atmospheric concepts are 
reviewed. The model and its computer program are described. 
The results of trend and sensitivity analysis reasonableness 
tests are presented. Finally suggestions for use are made. 
Many plots of mission life predictions are presented. The 
model computer program and sample input and output are also 


included. 
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eee NR ODUCT LON 


All current low-earth-orbit satellites are simply placed 
in orbit and then tracked as their orbits decay deeper into 
the atmosphere. This orbital decay is due to forces on the 
satellite produced by atmospheric drag and the  non- 
spherical earth. 

Modified exponential-type models for the atmosphere 
produce sufficiently accurate drag predictions to track the 
satellites. 

New concepts in scientific and military uses of space 
require greater accuracy in predicting atmospheric drag 
effects then produced by the modified exponential models. 
These new uses include the Aero-assisted Orbital Transfer 
Vehicle (AOTV), orbiting space stations, the Trans 
Atmospheric Vehicle (TAV), and the strategic defense 
initiative satellite constellation. 

The environment of all of these vehicles is called the 
thermosphere, an atmospheric layer studied in depth only 
recently. An accurate empirical model of the thermosphere 
has also been produced. 

One idea under investigation at the Naval Surface 
Weapons Center is the placement of many low-earth-orbital 
defensive satellites in precisely controlled and maintained 


orbits. This set of defensive satellites is termed a 


constellation. All low earth-orbital-satellites drift from 
their initial orbits due to perturbating accelerations 
caused primarily by the aspherical earth and atmospheric 
drag. In order to be functional, the orbits of the 
satellites in an SDI constellation must be exactly 
maintained. The drift of a satellite must be cancelled out 
by using an onboard propulsion system to produce forces to 
offset the perturbating accelerations. 

One model, proposed by Dr. Dan Parks of the Naval 
Surface Weapons Center, calculates the mission lifetime 
expected of a satellite doing intrack micro-thrusting to 
overcome atmospheric drag effects. The end of mission life 
corresponds to depletion of all the maneuvering fuel. This 
model is here-after referred to as the propellant longevity 
model. The propellant longevity model incorporates a Bessel 
function expansion for the fuel-mass decrement on each 
revolution. However, the model does not calculate the fuel 
necessary to cancel out orbital drift due to the aspherical 
earth. 

This thesis gives a verification of the low-altitude 
earth satellite propellant longevity prediction model 
developed by Dr. Parks. The theory for the model 
development itself is presented in Naval Surface Weapons 
Center Technical Report (NSWC TR) 83-243 (Parks [Ref. 1]). 
Computer code for the model was developed at NSWC under the 


direction of Dr. Parks. In addition to the computer code of 


the model, several other NSWC publications with applications 
to the model were forwarded to the Naval Postgraduate 
SCHOO]. 

The Thesis is organized as follows. Chapter I presents 
a general introduction to the propellant longevity model and 
discusses the changes made to the computer code for use at 
the Naval Postgraduate School. A basic description of the 
orbital mechanics necessary to operate the model is provided 
in Chapter II. Chapter II also presents an overview of the 
atmospheric physics necessary to understand the density and 
compositional variations in the thermospheric layer of the 
thermosphere. A description of the propellant longevity 
model is given in Chapter III. Instructions for operation of 
the computer code are also included there. Chapter IV 
contains the eee of the verification process sensitivity 
analysis and reasonableness tests. A summary is given in 
Chapter V. Appendix A contains the propellant longevity 
Fortran computer code for the IBM 3033. Appendix B shows 
input sets and their corresponding outputs. Graphical 
presentation of the tabular results is included. 

The coding of the model by NSWC was completed in Fortran 
for the Control Data Corporation (CDC) 6700 computer. 
Because of software and hardware differences, extensive 
modifications to the model coding were required to install 
the program as an operational model on the IBM 3033 computer 


at the Naval Postgraduate School. To optimize the utility 


of the model, some adjustments were also made to the 
original computer coding of the propellant longevity model. 
Great care was taken when modifying the code to preserve the 
model theory as contained in NSWC TR 83-243. 

The alterations to the computer program fall into two 
categories: changes to transport and implement the model on 
our IBM 3033, and changes made to mainstream the. model 
verification process. Changes to transport the model from 
the CDC environment to the IBM environment involved some 
logical revision and syntax, hardware word-size 
considerations. 

The hardware word-size problems _— extensive. The CDC 
6700 uses a 60-bit word size for single precision 
calculations, whereas IBM's standard is a 30-bit word. The 
ates STE longevity model requires 60 bit precision. Thus 
all real numbers, library functions, and input constants 
were promoted to double precision accuracy. Another 
precision-related problem involved the use of fractional 
exponentiation. All one-half exponents were changed to call 
the SQRT function to increase precision and speed. 

Syntactical changes were mostly related to the software 
environment transportation from CDC Fortran to IBM Fortran 
7 de Library functions calls, and input and output format 
statement changes, were the bulk of the syntactical changes 
made. (Examples of library function changes are the ARQTAN 


function changed to the ATAN2 function, the BESI Bessel 


function Bani changed to the MMBSIR Bessel function call and 
so forth). 

Examples of changes in the logic to in order to 
transport and implement the model code include using the 
Fortran 77 block IF-THEN-ELSE statement for the CDC Fortran 
IF statement and the use of DIMENSION statements in each 
subroutine that is processing a particular array. 

The second category of changes to the program addressed 
mainstreaming the model verification process. The model 
code as down-loaded from the NSWC tape contained little 
code-internal documentation. Looking towards future 
implementation and maintenance of the model code, extensive 
documentation was added as code internal comments. Most of 
this documentation was based on the research literature of 
NSWC. Extensive reference to the SUELSRE publications is 
mace in the code documentation itself. 

Two changes were made in the way the code solves the 
model in order to implement the theory with greater 
precision. The first change turned on the subroutine THRST 
which was designed to allow for preplanned changes to the 
semi-major axis of the orbit at the expense of onboard 
maneuvering fuel. (The subroutine code down-loaded from 
tape did not return the fuel consumed for the orbit adjust, 
record the resultant change in the semi-major axis, or 
receive the correct intended adjust parameter. All these 


functions are performed in the IBM 3033 version.) 


A second change in the code alters the way in which the 
satellite velocity relative to the atmosphere is computed. 
In the original code this parameter had to be entered as 
part of the input file and was assumed to be constant. 
Since this relative velocity is actually a function of other 
transient orbital parameters, recomputation is now 
automatically accomplished for each satellite revolution 
based on current values of the changing orbital parameters. 
In place of the satellite-to-atmospheric-relative-velocity 
input, the user now inputs the atmospheric-to-earth relative 
rotation rates. The physics associated with this change is 
discussed later_ on, in the section describing the 
atmospheric model. 

A few other changes to the propellant longevity code 
customize the output, or state the input element initial 
conditions. Several control inputs now allow for a range of 
model predictions. These are explained in Chapter III. 

The actual model verification is Seconpil shed in three 
steps. First does the model address the original problem? 
Second, does the model meet the. criteria of common sense? 
Finally, how well does the model solution compare when 
tested against real-world data? 

No real-world data sets exist for testing the propellant 
longevity model. The only low-earth-orbit satellite having 
a propulsion system to maintain orbital parameters is highly 


classified under the auspices of the United States Air 
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Force. The data from this program is not available in an 
unclassified state. In lieu of actual data, several 
additions to the program code were made to assist in the 
model verification process. The major adjustment is the 
addition of a feature to allow for a sensitivity analysis. 
The subroutine SNALYS, allows the user to linearly scan any 
combination of the 20 input elements, from a user specified 
minimum to any maximum, with user specified granularity. 
Ge Seutput of the Sensitivity analysis can be descriptive 
and tabular, or graphical in nature. The SNALYS subroutine 
allows for a plot of model prediction against the entire 
range of the input elements. Results of input element 
sensitivity analysis are presented in the Chapter IV. The 
sensitivity analysis feature allows for optimum orbital- 
element and satellite-ballistic-coefficient design. 

Thought was given to possible model installation on 
other computers. The model is large and requires about 20 
minutes of dedicated central processor time to complete an 
average run. As such the model is rather limited to 
mainframe computers. Initial calculations based = on 
processor time indicate an IBM AT microcomputer with the 
requisite 80287 co-processor and additional RAM memory might 
require about ten hours to complete an average run. 

This chapter presented an overview and motivation for 
the propellant longevity model. Changes made to the 


computer program at the Naval Postgraduate School were also 


id 


discussed. The next chapter presents the necessary 
background in orbital mechanics and atmospheric concepts to 


operate and understand the model. 


J 


II. BACKGROUND 


This chapter presents the orbital mechanics’ and 
atmospheric concepts required for the understanding and 
operation of the propellant longevity model. The orbital 
mechanics presented in the first section discusses only the 
orbital element sets, concepts, and coordinate systems 
specific to the model. In the second section of this 
chapter, Thermosphere specific atmospheric concepts’ and 


cycles are presented. 


A. ORBITAL MECHANICS 

This section covers the orbital mechanics for operating 
the mass decrement model. If the coordinate and orbital 
element systems discussed are familiar, this basic orbital 
mechanics section can be skipped. Orbital mechanics 
references are cited in the text. 

Early celestial mechanics models were based on 
increasingly accurate Sens: As observations and 
record keeping improved so did the empirical model fit to 
the observations. The model that fathered modern orbital 
mechanics was the empirical ences derived by Johannes 
Kepler, 1571-1630. In 1619 Kepler formulated three laws 
(models) of planetary motion. These models were based on 
life long observations recorded by Tycho Brahe, 1564-1601. 


Kepler's models were formulated as follows: 
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1. The orbit of each planet is an ellipse with the sun at 
one focus. 


2. The line joining a planet to the sun sweeps out equal 
areas in equal times. 


3. The square of the period of a planet is proportional 
to the cube of its mean distance from the sun. 


Later, Sir Isaac Newton, 1643-1727 validated Kepler's 
laws by Arsoeionliusiniens, their equivalence with his law of 
universal gravitational attraction. 

Newton's model for gravitational attraction states that 
two bodies attract with a force that is proportional to the 
product of their masses and inversely proportional to the 
square of the distance between them. Symbolically, the 


force on the mass m due to the mass M is given by: 
ig 
= (1) 


where G is the universal gravitational constant and r is the 
vector from the center of mass of M to the center of mass nm. 
: is the magnitude of the distance vector r. 

A body experiences an acceleration due to the resultant 
of all the forces acting on the body. If only a 
gravitational force acts on the body, then the acceleration 
of the body is given by: 

G(M+m) — 


rN A ee ee (2) 
Ls 
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In this model, r" is the acceleration vector on m caused by 
gravitational attraction of M on m. The remaining variables 
are defined as in the previous model. If M is much greater 
then m (as in an earth satellite system), G(M+tm) is very 
nearly GM. Another constant can then be defined called the 


gravitational potential parameter u: 
u = Gm (3) 


The gravitational acceleration of a satellite in earth 


orbit can then be expressed as: 
r? + u/r? *r = 0 (4) 


The gravitational field is conservative. This means 
that an energy constant of motion can be defined for a body 
moving in a gravitational field. For a two body, an- 
atmospheric, elliptic orbit system the specific mechanical 


energy is: 
F = vy?/2 - u/r (5) 
In this model ¢€ is the specific mechanical energy, v; is 


the satellite tangential speed, and u and fr are as 


previously defined. Note that E is the sum of the kinetic 
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energy per unit mass of the satellite and its potential 
energy per unit mass. 

The model, (5), requires satellites at specific radii to 
have specific tangential speeds. As r increases, v must 
decrease to maintain € constant. This orbital energy model 
is important when modelling the effects of atmospheric drag 
which is the major topic of the next section. 

Newton's inverse square law of gravitational force 
requires gravitational field orbits to be elliptical in the 
shape. The mathematics of the ellipse is important in the 
description of orbits and Figure i presents some of the 
important parameters. 

The parameters shown in Figure 1 are defined as follows: 

a: the semi-major axis. 
b: the semi-minor axis. 


e; the eccentricity where 


e = SQRT (1 - (b/a)?}. 


C: the ellipse center. 
O: the ellipse focus (earth center). 
p: the semi-latus rectum. 


r: the distance from a position on the ellipse 
the focus 0O. 


vz: the true anomaly. 
rp: the radius at perigee (in an earth system). 


ra: the radius at apogee (in an earth system). 
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Figure 1. The Elliptical Orbit 


the apogee (point farthest from the focus). 


(> 


P; the perigee (point closest ee the focus. 

No real orbits occur in just two dimensions. To describe 
real orbits a third dimension must be added. A suitable 
reference must also be chosen. One of the most convenient 
reference frames is the geocentric reference depicted in 


Figure 2. In this system a fixed nonrotating earth forms 
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the center of the celestial sphere. The earth's projected 


north pole, south pole, and equator become the north 


the south celestial pole (SCP), 
celestial equator respectively. 


eo 


Celestial Equator 


celestial pole (NCP), and 


SCP 


Figure 2. The Geocentric Celestial Sphere 


To account for the apparent motion of the celestial 


sphere aS viewed from the earth, a reference point on the 
celestial sphere must be chosen. Traditionally this has 
been the First Point of Aries, describing the constellation 


18 


where the vernal equinox occurred when the celestial 
reference was defined. The vernal equinox does not now 
occur in the constellation Aries. Tradition still uses the 
first point of Aries as the Celestial Meridian (CM) of 
reference. Figure 3 illustrates the vernal equinox 
reference. The vernal equinox is defined as the time when 


the sun is directly over the equator at noon in the spring. 


VERNAL 
EQUINOX 


© 


SUMMER 
SOLSTICE 






WINTER 


EQUINOX 


Figure 3. The Vernal Equinox 


dS) 


Figure 4 is a representation of an orbit in real space 
using the geocentric’ systen. The points where the 
projection of the satellites orbit on the celestial sphere 
crosses the celestial equator are called the nodes of the 
orbit. The ascending node (N) is where the projected orbit 
crosses the celestial equator heading for the north 
celestial pole. The descending node is where the projected 
orbit crosses the celestial equator heading for the south 
celestial pole. The angular measurement from the reference 
celestial meridian to the ascending node along the celestial 
equator is called the longitude of the ascending node, and 
usually given the symbol OQ. The angular measurement from 
the ascending node to the point of perigee (P) along the 
projected orbit is called the argument of perigee, 
symbolized by w. Both the longitude of the ascending node 
and the argument of perigee range from 0 to 360 degrees. 
The angle from the plane of the celestial equator to the 
mene of the projected orbit is called the inclination and 
denoted by i. Values for the inclination range between 0 
and 180 degrees. 

Inclinations of 0 and 180 degrees are termed equatorial. 
In equatorial orbits the longitude of the ascending node is 
not defined. Orbits with an inclination of 90 degrees are 
called polar. Prograde orbits are those with an inclination 
of between O and 90 degrees and retrograde orbits have an 


inclination range of 90 to 180 degrees. 
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NCP 





Figure 4. A Geocentric Orbital System 


Historically, it was difficult to tell much about the 
orbit specific parameters from the earth vantage point. 
However, a value that could be measured was the time for one 
complete orbit. If the radius vector moved uniformly along 
the projected orbit with the average rate of mean motion n 
(n is also referred to as the mean angular velocity) at an 


epoch T, it would be at a position M called the mean 
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anomaly. The mean anomaly has a range of 0 to 360 degrees. 
It is measured from the perigee along the projected orbit. 
The mean anomaly corresponds to positions on the orbital 
ellipse. An eccentric anomaly (E), is defined as the angle 
between the projected position and the point of perigee 
subtending at the center of the ellipse. Figure 5 


illustrates eccentric anomaly. 





Figure 5. Orbital Anomalies 


oe 


The difference between the eccentric anomaly and the mean. 
anomaly is important. The mean anomaly is a mathematical 
concept corresponding to an imaginary position of a 
satellite past the perigee if it moved at an average angular 
velocity. The eccentric anomaly is the actual projected 
angular position of a real satellite on the celestial 
sphere. The eccentric and mean anomalies are Scan in 
circular orbits (e = 0). The eccentric and mean anomalies 


are related by Kepler's equation: 
M = E - e*sin E (6) 


The solution to model (6) has been the source of much 
investigation. By a solution it is meant for a given M, 
determine e and E. In 1817 Friedrich Wilhelm Bessel, 1784- 
1846 invented Bessel functions to provide an infinite series 
approximation for the solution of the model (6). Newton 
developed a converging algorithm for its solution. The mass 
decrement submodel of the propellant longevity model uses 
Newton's method for the solution of Kepler's equation. 

All anomalies are concerned with the measunenene of the 
pesteton Cf a body in orbit. All anomalies are measured 
angularly from the point of perigee in the orbital plane, 
and vanish at perigee. 

When e and E have been found from Kepler's equation (6), 


then using the equations for an ellipse in Taff [Ref. 2:p. 
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67], the remaining elliptical parameters can be solved for 
when any one of a, b, rp, or ra is known. 

The minimum number of parameters needed to describe an 
orbit is called an orbital element set. Five parameters are 
needed to describe an orbit. One more is needed to position 
an object in the orbit. Many orbital element sets have been 
used and are in use. The classical orbital elements are: 

semi-major axis (a) 

eccentricity (e) 

inclination (i) 

longitude of the ascending node (2) 

argument of perigee (wy) 

time of oes ee eaee or epoch (T). 
The Keplerian set replaces the time of perigee passage with 
the mean anomaly (M). The orbital element set used in the 
mass decrement submodel is the Kozai mean element set, very 
Similar to the Keplerian set. The name "Kozai mean element 
set" refers to the Kozai perturbation approximations used to 
solve for the Keplerian orbital element’ set. The 
perturbation theory approximation of Y¥. Kozai is discussed 
later. 

Up to now all orbital considerations have been limited 
to the classical two-body problem. The two body problem 
limits the forces on the orbiting body to the gravitational 
force. The two-body problem also requires all masses be 


concentrated at a single point. A real earth satellite 
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experiences many forces resulting in many discrete 
accelerations. The theory of perturbations was developed to 
account for the motion of real satellites. 


In perturbation theory the original orbit, with only the 


point-mass gravitational acceleration, is called the 
osculating orbit. Orbital elements describing the 
osculating orbit are the osculating element set. The 


osculating orbit is then adjusted periodically to account 
for the action of all the nonpoint-mass gravitational 
accelerations. 

Accelerations result from forces. Examples ae 
perturbating forces on an earth satellite are: oblate earth 
gravitational, atmospheric drag, sun gravitational, moon _ 
gravitational, solar radiation pressure, and satellite motor 
eEnNruUSt. 

To remain ila a closed eels pay the point-mass 
gravitational force must be the dominant force on the 
satellite. For low earth orbit satellites, the oblate earth 
gravitational force is the next greatest force and is termed 
the dominant perturbation. Satellites below 600 kilometers 
in altitude experience an atmospheric drag force of the same 
order of magnitude as the oblate earth force. The 
atmospheric drag and oblate earth effects are so large for 
low earth orbit satellites that the effects of the other 
perturbations (except motor thrust) are .lost in the 


background error. 
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Central to the problem of perturbation theory is the 
that the three body problem has never been solved 
mathematically in general closed form. Most modern 
approximations involve power series or Legendre polynomial 
expansions. The state of the art in orbit perturbation 
prediction only uses only the first-order terms of these 
expansions, termed first-order effects. An exception is a 
theory advanced in 1959 by Y. Kozai. His approximations are 
computationally efficient and allow for the use of more of 
the oblate earth effects (called J2, J3, and J4 terms and 
discussed later). Kozai's theory is flawed in the higher- 
order terms due to the introduction of mathematical 
singularities. Nevertheless, proponents for its use argue 
that it does work and predicts the oblate earth effects toa 
reasonable degree of accuracy. The mass decrement submodel 
of the prepa iene longevity model presented in this thesis 
uses Kozai's theory. 

The gravitational potential parameter u was presented as 
Singularly determining the gravitational potential in the 
earlier development of the gravitational acceleration model. 
In earth systems, the gravitational potential is a complex 
function of the earth's mass distribution. 

The earth is not a point mass. Neither is it a 
homogenous sphere. It's true mass distribution is complex. 


The real distribution of matter on the earth forms a body 
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termed the  geoid. This geoid determines exactly the 
gravitational potential function. 

The determination of the geoid is based on ground and 
satellite measurements. Many military programs depend on 
the precise position of a point on the surface of the earth. 
The military has led the way in geoid determination in the 
past two decades. 

Geoidal studies reveal important facts about’ the 
aspherical shape of the earth. A fluid spinning sphere 
should exhibit a bulge at the equator and flattening at the 
poles. This is termed oblateness. The oblate deviation is 
the dominant deviation for the earth as ae_=e sphere. 
Measurements of the earth's oblateness indicate that it is 
more oblate than the current spin rate could produce by 
itself. This excess oblateness Pancakes two factors: the 
earth's core is more solid than previously thought, and the 
earth's spin has slowed. 

The earth deviates in other ways from being a sphere. 
Since all deviations can be expressed in terms of sines and 
cosines they are referred to as spherical harmonics. 
Harmonics are divided into categories. Zonal Havkontes are 


symmetric about the earth's spin axis and are not longitude 


dependent. Sectorial harmonics are only longitude 
dependent. Harmonics dependent on both longitude and 
latitude are termed tesseral. Odd numbered harmonics are 
antisymmetric about the equatorial plane; even numbered 
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harmonics are symmetric about the equatorial plane. 
Harmonics are also classified by number. The number refers 
to the subscript of the J coefficient (defined later). 

Bate et al., [Ref. 3], offers a simplified development 
of the gravitational potential function. If the 
gravitational potential function 6 of a geoid is known, the 
acceleration due to this gravitational function is modeled 


bys 
rr" = Vo (7) 


One common expression for the gravitational potential or 


geoid function @ is: 
Uu i Te Nl (8) 
eh 2, J. (—=)eebaesomele) 


In this model u is the gravitational potential parameter, Jn 
is an experimentally determined J coefficient, r, is the 
radius of the earth, Pn is a Legendre polynomial of order n, 
L is the geocentric latitude, and r is the radius of the 
satellites position. Taking the gradient of ante potential 
function yields an expression for the acceleration due to 
the gravitational potential. Bate et al., [Ref. 3:pp. 419- 
423] details the resultant potential induced acceleration to 


the first seven terms. 
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Recent studies have determined the geoid function out 
through the J44 term. Values used by the mass decrement 
model by the subroutine GEOP are 1082.76*1076, -2.55*1076, 
and 1.56*1076, for the J2, J3, and J4 terms respectively. 
The J2 and J4 terms are even numbered, and refer to even 
harmonics symmetric about the equator. The J3 term is odd 
numbered, and refers to odd harmonics which are 
antisymmetric about the equator. The J2, J3, and J4 terms 
are all zonal. 

Current perturbation theory expresses @ - to the J2 
term. The theory advanced by Kozai in 1959 appears to use 
the J2, see and J4 coefficients in a first-order 
perturbation theory approximation (because singularities 
occur in terms higher-order than J2). Taff [Ref. 2:pp. 332- 
342] takes strong exception to the theory. 

The subroutine GEOP in Se propellant longevity model 
uses Kozai's theory to calculate the geoidal effects on 
orbital parameters through the first three zonal 
coefficients (J2, 33, and J4 terms). Their effect is 
Calculated once each orbital revolution, and the Kozai mean 
element set is adjusted accordingly. While using the model, 
it is important to understand that as the altitude of the 
satellite decreases, The aspherical earth perturbations | 
become much greater. The orbital element drift due to the 
aspherical earth is thus much faster in lower altitude 


orbits. 
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In summary, the computer program of the propellant 
longevity model uses three orbital element sets in a 
geocentric reference. The input orbital element set is the 
Brouwer mean element set of the Naval Space Surveillance 
System (NAVSPASUR)1! system. Chapter III and Appendix A 
explain the input elements in detail in the input 
description. The Brouwer mean elements are first converted 
into a Brouwer osculating set. These osculating orbital 
elements are then transformed into a Keplerian-like Kozai 
mean set. Half of the computer program of the propellant 
longevity model transforms the Brouwer mean element set into 
the Kozai mean set. These element set transformations are 
used to prepare real-world orbital element sets for use in 
the propellant longevity predictions. Once the computer 
program has solved for the Kozai mean element set all 
perturbation effects geste in adjustments to the Kozai mean 


set. Adjustments are made once per revolution. 


B. THE ATMOSPHERE MODEL 
Central to the validation of any model iss-~— an 
understanding of the system to be modeled. The model under 
investigation determines the propellant longevity for low 
earth orbital satellites doing in track micro-thrusting to 
ltfhe NAVSPASUR system publishes a semiannual report on 
the 1000+ orbiting objects they track in The Satellite 
Situation Summary. The orbits are described ina five card, 
80 column, Fortran style input format. The report is 


avallable on request from the Commanding Officer, Naval 
Space Surveillance System, Dahlgren, Virginia. 
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overcome the effects of atmospheric drag. The fuel required 
by a particular rocket motor on a satellite to perform a 
known amount of work can be computed from the rocket 
equation. The problem occurs in computing the amount of 
work required. 

To compute the required work by the satellite propulsion 
system to overcome atmospheric drag requires knowledge of 
the atmospheric environment of the satellite. For low earth 
satellites the atmospheric environment is the thermosphere. 
The physics of the thermosphere is discussed in this 
section. Before discussing the thermosphere, some 
preliminary thermospheric properties are presented. 

For low earth orbital satellites, the drag-induced 
change in the orbital elements is significant. The oblate 
earth (referred to as the "J2 term") orbital perturbation is 
normally the dominant perturbation term. The perturbations 
due to drag are of the same order of magnitude as those due 
to the oblate earth for orbits below 600 kilometers. 

The drag force caused by the satellite moving through 
the atmosphere always acts against the motion of the 
satellite. Neglecting the lift effects due to the 
satellite's shape and attitude, the drag force Patces work 
to be done on the atmosphere by the satellite. 

From Bate e al., [Ref. 3:pp. 15-16] a satellite in orbit 


has a constant specific mechanical energy §& given by 


Spit 


E = = -¢ (9) 


where € is the mechanical energy per unit mass, vy the 
satellite speed, u the gravitational potential constant, and 
r the radius of the orbit. 

The above model says that the mechanical energy per unit 
Mass of a body in an an-atmospheric, two body, circular 
orbit is the sum of its kinetic energy per unit mass and 
gravitational potential energy per unit mass. 

Conservation of energy requires that the work done on 
the atmosphere by the satellite must reduce the orbital 
energy of the satellite. The differential element of work 
done on the atmosphere by the satellite results in a 
differential reduction in the specific kinetic energy of the 
satellite. Assuming the specific mechanical energy is a 
constant, the radius, r, of the orbit must therefore 
decrease. As the radius decreases, the velocity tends to 
increase in order to maintain the an-atmospheric 
conservation of specific mechanical energy. 

Abe non-circular (real elliptical) SEDs? the 
atmospheric effects are greatest at perigee (Figure 1). 
Orbital mechanics requires a differential reduction in 
orbital velocity at perigee to result in a reduction in the 
semi-major axis and orbital eccentricity. The major effect 
of transatmospheric satellite motion is to reduce the semi- 


Major axis and the eccentricity of the elliptical orbit. 
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Taff (Ref. 2:pp. 352-353] develops the mathematical 
equations for the drag-induced changes in the semi-major 
axis, eccentricity, argument of perigee (Figure 4), and mean 
anomaly (Figure 5). The changes in the argument of perigee 
and the mean anomaly are shown to be periodic. The net 
effect on the orbit due to atmospheric drag is to 
circularize the orbit closer to the earth, with little 
change to the other orbital elements. 

An accepted model for the magnitude of the drag force as 
presented in Taff [Ref. 2:p. 352] is: 


ae SS 2 


Fn = — V 


2m 1) 


In this equation, Fp is the drag force magnitude, A the 
effective satellite cross-sectional. area, Cp the drag 
coefficient, m the satellite mass, p the atmospheric 
density, and Vre, the velocity of the satellite relative to 
the atmosphere. 

The upper atmosphere is dynamic. The changes that must 
be modeled to achieve a valid prediction for the drag force 
are those that effect the drag force factors. The factors 
in the drag force model that are dependent on the atmosphere 
are density, a drag coefficient (dependent on satellite 
physical ‘and fluid-flow characteristics), and the relative 


velocity of the satellite in the atmosphere. 
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To model the atmospheric effects on density, drag 
coefficient, and relative velocity, some background in upper 
atmosphere physics is needed. The descriptive presentation 
here is carefully developed to be relative to the above 
three drag force dependent factors. A more comprehensive 
treatment of atmospheric physics can be found in a work of 
Fleagle [Ref. 4]. An axeerTene technical treatment of 
thermospheric dynamics can be found in the works of Jacchia 
(Refs. 5,6,7]. 

Jacchia formulated the most widely used thermospheric 
model currently employed by professional orbital analysts 
(known as the Jacchia J60 model or simply J60). His latest 
comprehensive empirical thermospheric model (Jacchia J77 or 
J77) is many orbital analysts' preferred refinement for 
modeling the atmosphere. 

The Jacchia J77 =model is highly tabular, uses 
considerable computer time, and requires an extensive input 
file. However, these negative aspects are balanced by the 
orders-of-magnitude increase De) accuracy of density 
predictions. Jacchia [Ref. 7:p. 2] presents a table of 
residuals for predicted versus observed densities for 
various satellites, as well as time intervals for the 
Jacchia J77 atmosphere model. All density residuals are 
much smaller than those produced by exponential models in 
the same time intervals. When extremely accurate density 


predictions are required in some submodel, the use of the 
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Jacchia J77 atmosphere model may well be worth the 
additional computer time and expense. 

A simple model of the atmosphere considers it as a 
homogeneous mass of constant density. A refinement is made 
by modeling the density as decreasing exponentially with 
altitude. The exponential model predicts pressure and 
density will fall one tenth for every ten kilometer increase 
in altitude. Prior to 1960 this exponential model was used 
extensively for upper air. calculations involving meteor 
trails. With the advent of low earth orbital satellites, a 
greater degree of accuracy in predicting atmospheric effects 
is required. The predicted satellite decays from the 
exponential model are very different from the observations. 
The Jacchia J60 atmosphere model was formulated to improve 
decay predictions. It was based on measurements form high 
altitude balloons, sounding rockets, and ground stations. 

The atmosphere is divided roughly into concentric 
shells, each with specific properties. The properties used 
in this subdivision are after Fleagle. [Ref. 4] 

From the early days of atmospheric studies, temperature 
and moisture have been recorded. The division system 
described by Fleagle [Ref. 4] is based primarily on a height 
versus temperature profile. 

The height versus temperature profile is graphed in 
Figure 6. The temperatures and altitudes are those of the 


1975 U.S. standard atmosphere. This temperature behavior 
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Figure 6. Atmospheric Temperature Profile 


approximates general conditions (large transient variations 
have been observed). Figure 6 also shows the boundaries 
between concentric shells, called pauses. The following 
description of the atmospheric layers is presented in 
Fleagle [Ref. 4:pp: 79-84]. 

The region of the atmosphere closest to the earth is the 


troposphere and contains eighty percent of the atmospheric 
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mass. This region is characterized by near’ linearly 
decreasing temperature with altitude. The dynamics are due 
to earth surface energy transitions. All of the earth's 
weather occurs in this’7 region. “The height of the 
troposphere is variable and depends on both latitude and 
season. Along the equator the maximum altitude is 16 to 18 
kilometers. At the poles the troposphere is highly 
seasonal, being nearly absent during the winter months and 
extending to eight to ten kilometers in the summer. Capping 
the troposphere is a region of constant temperature with 
altitude, called the tropopause. 

Next in the series of atmospheric layers is the 
stratosphere. The stratosphere is characterized by linearly 
increasing temperature with height. The physical dynamics 
are due mainly to the absorption of solar ultra-violet 
radiation by the ozone. The maximum altitude of the 
stratosphere is about 50 kilometers where it is capped by an 
area of local maximum temperature termed the stratopause. 

Above the stratopause is the mesosphere. This 
atmospheric layer exhibits linearly decreasing temperature 
with height. The energy physics of the mesosphere are not 
well understood. Bulk mixing of the gasses still dominates 
diffusion as thermodynamic processes. Any molecular 
dissociation of atmospheric molecules produced by radiation 


absorption is rapidly consumed by chemical recombination of 
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these dissociated species. The mesosphere is topped by the 
turbopause at an altitude of about 90 kilometers. 

Above the turbopause the gross physics of the atmosphere 
changes drastically. MThis equisn is the environment of all 
current low earth orbit satellites. Termed the 
thermosphere, it extends from the turbopause to an altitude 
of 1000 to 2500 kilometers. Here bulk molecular mixing is 
dominated by diffusion. Lighter gasses diffuse upwards and 
heavier ones remain at the lower altitudes. Disassociation 
of molecular species becomes increasingly important, this 
results in an increase of individual atomic’ species. 
Standard flow equations (Navier-Stokes) are extremely 
difficult to apply, and physical descriptions are more 
discrete (quantum mechanical). Flight by flow-supplied lift 
is not possible. The standard methods of computing drag 
coefficients do not apply. The ideal gas approximation is 
also not valid (P*V = n*R*T does not apply). Above the 
turbopause, density is an increasing function of 
temperature. Only about one millionth of the atmosphere 
remains above the tropopause. 

The stratification of gaseous components’ in the 
thermosphere is significant. The drag coefficient has been 
shown to be composition dependent. Thus, as the gaseous 
composition changes the drag coefficient must be adjusted to 
reflect the observed drag-induced orbital decays (Jacchia 


(Ref. 7]). Figure 7 shows the dominant atmospheric species 
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Figure 7. Atmospheric Composition 


by altitude with the accepted drag coefficients. The drag 
coefficients are predictions obtained analyzing the drag 
coefficient needed to produce the observed drag effects on 
orbital decay for a satellite at a particular altitude. 
Thus the drag coefficient is a proportionality constant in 


the drag force model. 
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For altitudes below 110 kilometers, molecular nitrogen 
is the main atmospheric gas. Above this region the 
dissociation of the molecular oxygen by ultra-violet light 
absorption causes atomic oxygen to become increasingly the 
main gaseous constituent. At an altitude of 500 kilometers, 
helium replaces atomic oxygen as the main atmospheric gas. 
At altitudes above 1500 kilometers atomic hydrogen becomes 
the dominant gas (see Figure 7). 

Traditional research has assigned a drag coefficient of 
2.20 to altitudes between 200 and 400 kilometers. As the 
altitude of the satellite increases above 500 kilometers, 
corresponding to satellite movement through helium, a drag 
coefficient of at least 3.00 must be assigned to maintain 
agreement with observation (see Figure 7). 

Jacchia [Ref. 7] lists the observed variational changes 
in density and composition in ‘the thermosphere. They are 
the two component solar activity related, semiannual, 
diurnal, seasonal latitudinal, and geomagnetic. 

The two component solar thermospheric variations are 
consequences of changes in solar ultra-violet (EUV) flux. 
The EUV flux changes periodically due to two solar cycles. 
These are the 27 day rotational or disk cycle and the 11 
year sunspot or active area cycle. 

EUV flux is not measurable from the ground due to almost 
complete absorption in the atmosphere. Empirical evidence 


has suggested a strong parallel of the EUV flux and the 
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ground measurable 10.7 centimeter flux (F10.7) which has 
been used in some successful thermospheric models. The 
Jacchia 1960 model used the 10.7 centimeter flux as an input 
to model thermospheric density resulting in an improvement 
of decay predictions over the exponential model. (Indeed, 
it is only recently that there has been a requirement for 
greater accuracy in drag predictions.) 

The thermosphere responds to increased EUV flux by 
absorbing more energy, causing increased heating. This 
-increased heating causes an increase in temperature and 
density. The thermospheric response to changes in the EUV 
flux is nonlinear. Increases in EUV flux due to disk 
(rotational) effects produce more thermospheric heating then 
a like increase in flux due to active area effects. 
However, Tee noticed in 1973 that the disk component of 
the EUV flux can be linearly related to an average 10.7 
centimeter flux. He recommends using a 71 day average 
corresponding to three solar rotations [Ref. 7]. Tables for 
the 10.7 centimeter flux have been developed. ? The Jacchia 
1977 atmosphere incorporates the average 10.7 centineéter 
flux as an input, used as an indicator of the thermospheric 


heating capacity of the sun's radiation. 


2Orr [Ref. 8:Appendix B] lists values compiled and 
predicted by the NASA Marshall Space Flight Center. Values 


range from a minimum of 73.3 to a maximum of 242.9. The 
units fos the LOe 7. centimeter flux are 10722 
watts/meter/Hz. 
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The semiannual temperature and density variation has 
been modeled in many ways. Jacchia [Ref. 7:pp. 47-50] fits 
a modified harmonic empirical model to the biannual, 
Sinusoidal cycle of observed density. Maxima occur in April 
and October. The Jacchia 1971 and 1977 models for the 
semiannual density variation are based on observations of 12 
years of satellite data. The observed density variations 
are thought to be a consequence of the earth-sun 
orientation. The mechanism for the observed effect is 
unknown, but may be related to the maximum equatorial 
heating at the equinoxes. Equinoxes (sun at maximum 
altitude over the earth's equator) occur in late March and 
September. Figure 8 illustrates the semiannual bulge in 
temperature and density observed just after the equinoxes. 
Due to the earth's oblation, the bulk mass of the atmosphere 
is concentrated at the equator. The most direct exposure of 
the equatorial atmosphere occurs at the equinoxes. The 
temperature and density thus peak shortly after periods of 
maximum heating. 

Seasonal-latitudinal thermospheric variations are again 
earth-sun orientation related. Variations in density and 
composition are dependent on the season and latitude. The 
Jacchia 1977 model fits modified harmonic empirical models 
to observed data. The seasonal latitudinal variations 
interact strongly with the semiannual variation, complicat- 


ing the modeling. The winter helium bulge is the best known 
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Figure 8. The Semiannual Atmospheric Density Bulge 


effect of this variational factor. Figure 9 shows the 
general shape change of the temperature and pene bulge 
with season. Note that during the June solstice the bulge 
is centered at about 45 degrees north latitude whereas 
during the December solstice the bulge is centered at about 


45 degrees south latitude. 
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Figure 9. The Seasonal-Latitudinal Atmospheric 
Density Bulge 


The diurnal bulge is shown in Figure 10. It is caused 
by the daytime heating of the side of the atmosphere facing 
the sun. Lagging the actual (most direct) solar exposure by 
about four hours, the heated side shows a marked increase in 
temperature and density. Many current models track and 


predict the position of the diurnal bulge. 
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Figure 10. The Diurnal Atmospheric Density Bulge 


Geomagnetic heating. of the thermosphere is not well 
understood. At times of high geomagnetic activity 
(corresponding to an index of 7) large heating eerents are 
produced. The combined effects of geomagnetic and EUV 
heating can produce changes in density of three orders of 


magnitude at 600 kilometers (Jacchia (Ref. Sip. 82]}). 
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In summary, the thermosphere is a very dynamic region of 
the atmosphere. Many of the observed changes require 
further research to complete analytical models. The current 
modified exponential empirical models can be orders of 
magnitude in error in their density predictions. At higher 
earth orbits orders of magnitude error in density is of 
little consequence. 3 

The model chosen for density and compositional effects 
in the atmosphere depends on the importance of the short- 
term cyclic effects. Computations for orbits of months in 
Guration and less than 600 kilometers in altitude will be 
subject to large errors using the simple empirical models 
(Jacchia 1966). For these short missions the more accurate 
Jacchia 1977 model seems more appropriate. (Performance of 
the Jacchia 1977 empirical model has been remarkable with 
density errors rarely exceeding two percent.) Missions 
lasting for decades may be able to use the computationally 
more efficient Jacchia J60 model to predict very long term 
drag effects. 

The one drag force dependent factor not yet discussed is 
the satellite-to-atmosphere relative velocity. This 


relative motion is influenced by three gross factors: local 


3Above 800 kilometers the drag term in the orbital 
decay model used by the Satellite Control Facility in 
Sunnyvale, California is very small. Its main function is a 
catch-all for errors in other perturbating terms. 
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winds, satellite tangential velocity, and aeriospiemc 
rotation. 

Thermospheric wind has not been measured. Mesospheric 
wind has been measured based on the observation on meteor 
tra@als. Fleagle [Ref. 4:p. 82] reports meospheric winds 
exceeding 0.150 kilometers per second. This may sound 
large, but must be compared with the satellite's velocity 
(which is on the order of 10 kilometers per second). 

Satellite tangential velocity is a major factor in the 
satellite-to-atmosphere relative velocity. When combined 
with the rotating atmosphere it accounts for the bulk of 
relative motion. 

Most earth orbit calculations simplify computation of 
perturbation-induced changes in the orbital elements by 
assuming a fixed earth. When considering drag a non- 
rotating earth assumption can lead to a severe error. The 
earth is not fixed. The atmosphere is pulled along with its 
rotation, and lags the earth rotation by a small fraction. 
Figure 11 shows the rotating atmosphere effect. For low 
inclination (prograde) orbits the atmosphere is actually 
following the satellite. In high inclination (retrograde) 
orbits the atmosphere is moving against the satellite. When 
‘comparing a typical 16 revolution per day satellite 
tangential velocity in a prograde atmosphere to the same 
satellite in a retrograde atmosphere, a 15 percent variation 


in relative velocity can occur. This is compared to the 
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Figure 11. Orbital Inclinations ina 
Rotating Atmosphere 


unknown (but probably less than two percent effect) 


atmospheric winds. 


of 


A model used by Parks [Ref. l:p. 4] to account for the 


rotating atmosphere separates the relative velocity of the 


satellite in relation to atmosphere into two parts. 


model is given by: 
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The 





“ee o= rE CVO (i) 


Here Vre) is the satellite to atmosphere relative velocity, 
vy is the satellite tangential velocity, and § is a factor 
accounting for orbital orientation (inclination) with 
respect to a rotating atmosphere (see Figure 11). 

Equation 11 in Parks [Ref. l:p. 4] models’ the 

inclination factor in the following way: 

r_ 2 
ele ieiomcos! i) (12) 

In this model, 6 is the rotation factor, lp the radius at 
perigee, Vp the satellite velocity at perigee, we the 
earth's angular rotation rate (assumed constant), i the 
orbital inclination, and A the ratio of rates of rotation of 
the atmosphere to — earth's rotation. 

The model for § used in the propellant longevity model 
neglects the effect of the atmospheric wind, but this should 
introduce only a small error. Recent atmospheric work, 
Jacchia [Ref. 7], suggests a significant gravity wave 
induced motion in the upper atmosphere. As more experimen- 
tal information on this bulk motion becomes available, the 
incorporation of the gravity wave wind may become 
increasingly important. 

The parameter A used by the propellant longevity model 


in the satellite velocity multiplier submodel is probably 
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not constant. Atmospheric theory suggests it may be 
dependent on both altitude and latitude. The range of 
values for the A parameter is assumed to be from a maximum 
of about 1.0 (at low altitudes) to a minimum of around 0.9 
(at high altitudes). 

The coding of the theory in NSWC TR 83-243 for satellite 
fuel mass decrement did not compute <6 internally. It was 
user supplied input in the original version. As part of the 
mass decrement model validation at the Naval Postgraduate 
School, code internal computation of 6 wae added. This 
allows 6 to automatically adjust the satellite-to-atmosphere 
relative velocity Vye, for the results predicted by the 
theory as perturbation-induced changes to the inclination 
OCCU. For a long mission life that exhibits significant 
inclination drift, the internal computation of should add 
accuracy to cine mass decrement submodel. In place of the 
original user input for 6, the Ais now specified by an 
input file. The input file default value for A is 0.9. 

A model atmosphere used for drag computations must be 
able to correctly predict the effects on density, drag 
coefficient, and satellite-to-atmospheric velocity. There 
are many variations in the atmosphere that affect density. 
The sum of all these effects creates a density bulge in the 
atmosphere, the location and shape of which follow maximum 


heating. 
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The drag coefficient is atmosphere composition 
dependent. Thus the model atmosphere must be able to 
predict the dominant species in the satellite environment. 

Satellite-to-atmosphere relative velocity is a function 
of satellite tangential speed, atmospheric rotation, and 
orbital inclination. 

The Jacchia J60 model atmosphere used in the propellant 
longevity model only includes the effects of diurnal heating 
on the density bulge. There are no computations of 
atmospheric composition. The newer Jacchia J77 model 
incorporates all the known density variations currently 
observed. The Jacchia J77 atmosphere also _ predicts 
thermospheric composition. 

The background concepts for the propellant longevity 
model understanding and operation have a presented. The 


next chapter describes the model and its computer program. 
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III. PROPELLANT LONGEVITY MODEL DESCRIPTION 


This chapter describes the propellant longevity model. 
First, a brief summary of the model theory is presented. The 
second part of this chapter describes the coding of the 
model as installed at the Naval Postgraduate School. 
Finally, the third section gives the model code compiling 
and handling instructions specific to the system at the 


Naval Postgraduate School. 


A. PROPELLANT LONGEVITY MODEL THEORY 

A detailed development of the propellant longevity model 
is presented by Dr. Parks [Ref. 1]. The propellant ° 
longevity model hinges on the prediction of the fuel mass 
used by a satellite to compensate for atncenhene drag. The 
mass decrement submodel predicts the amount of fuel consumed 
after each revolution of the satellite. 

The mass decrement submodel for a low-earth-orbit 
satellite in ang obiE Ee diurnal atmosphere is presented in 
Reference 1 as Equation 56. The diurnal density bulge is 
the only density and compositional atmospheric variation 
considered in the model. | 

The integrals of the mass decrement model (Equation 56) 
are next expressed in terms of Bessel functions of the first 


kind and order n as defined by: 
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271 
ieee) = = ! cos(nE) exp(fae cosE) dE (13) 


Here I, is the Bessel function of the first kind and order 
n, 8 the inverse density scale height (assumed constant), a 
the semi-major orbital axis, e the orbital eccentricity, n = 
0,1,2, ..-, 6, and E the orbital eccentric anomaly. All of 
these orbital mechanics terms were presented in Chapter 
II.A. The integral is to be evaluated over one revolution. 
The Bessel function expansion of the mass decrement 
model is given in equation 58 of Reference 1. This 
particular expansion is central to the model. Provided 
accurate values for other multipliers are produced by the 
model atmosphere, (accurate values for atmospheric density, 
satellite drag coefficient, and satellite-to-atmosphere 
relative velocity are essential) greatly increased accura- 
cies should result in predicting the fuel consumed to 
compensate for atmospheric drag. The main draw back to the 
use of the this type of Bessel function is the known 


Slowness of its convergence. 


B. COMPUTER CODE DESCRIPTION 

A detailed description of the model code and its input 
elements is contained in the code and included in Appendix 
A. An operational description is presented here. 

The computer code calculates the propellant mass 


remaining together with the cumulative mission life after 
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each revolution of the aoa When propellant fuel is 
Beeueted the computation terminates. 

Computer output is available in several forms and is 
specified by the user through the use of some control 
parameters. The user can select one of two modes; standard 
mode or a sensitivity analysis mode. The input elements used 
as control parameters for output and mode control are: 
maximum number of revolutions ITER1, sensitivity analysis 
granularity PSIZE, output print frequency ITER2, mode 
control PSENS, anomalistic period print control IET, number 
of scheduled orbit adjusts NOA, scheduled orbit adjust array 
IRV(K), and the amount of each adjust array DA(K). These 
terms are discussed below. 

The most basic form of output available is a tabulation 
of BrOpeTene remaining and accumulated mission life for 20 
specified Sea eee eal initial conditions. The 20 
initial conditions are specified in an input file. 
Accumulated mission life is expressed in modified Julian 
days (MID) 4 and revolution number. Frequency of output 
printing is selectable from a maximum of once per revolution 
to a minimum of once per 10000 revolutions by user 
specification of ITER2. 

A cap on runaway computer time is available by selection 
of ITER1. Satellites with small cross-sectional areas, high 

4x description of the modified Julian day time 


measurement is in Taff. [Ref. 2:p. 103]. The modified Julian 
day is reckoned from o2 Nov. 17 1858. 
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altitudes, and large amounts of maneuvering fuel can have 
very long orbital life-times. By defining a maximum number 
of revolutions to be considered in the calculations through 
the reasonable selection of ITER1, unnecessary expenditure 
of computer time can be avoided. 

As part of the model verification a sensitivity analysis 
option was added to the progran. It can be selected by 
setting PSENS = 1. The standard operating mode is selected 
if PSENS = O. In the sensitivity analysis mode, the 
subroutine SNALYS creates an array of elements selected by 
the user to be linearly scanned from a specified minimum to 
maximum with a given granularity. 

Granularity is specified by the selection of the control 
parameter PSIZE. A granularity of 1/100 (code specified as 
PSIZE = 100) will cause 100 executions of the propellant 
longevity code to complete fuel exhaustion. This 
specification can lead to excessive computer execution times 
for high altitude satellites. A granularity of 1/5 will 
cause five executions of the .code to fuel exhaustion. 
However, only five iterations from a given minimum to a 
given maximum may be too coarse to reveal the true trend of 
the satellite mission life. 

Each time the code is sequentially executed in the 
sensitivity analysis mode, the elements selected by the user 
for analysis are incremented linearly. The output in the 


sensitivity analysis mode is a table of mission lifetimes, 


a2 


resulting from the initial values of input elements not 
selected for analysis, and the current value of the elements 
selected for analysis. 

The output always echoes the initial input element 
conditions as read from a user supplied input file. Further 
choice of output and mode control parameters allow the user 
to tailor the output. 

Three more input control elements allow for the 
specification by revolution number of up to ten adjusts to 
the semi-major axis. The element IRV(K) is an array whose 
elements specify the orbital revolution number for axis 
adjust. The element DA(K) specifies the amount of change in 
the semi-major axis in kilometers for each scheduled 
adjustment. Finally, NOA specifies the total number of 
orbital adjusts to be performed. If the user does not | 
desire any changes to the semi-major axis, NOA must be set 
to O and no values for IRV(K) and DA(K) input to the file 
stream. 

An additional twenty. input elements must be specified by 
the weer in the input file. Five of these are satellite 
physical parameters, five are atmosphere specific 
parameters, and ten are the Brouwer mean orbital element set 
of the NAVSPASUR systen. 

The five satellite physical parameters are: propulsion 
motor specific impulse (SPI in seconds), initial mass of the 


satellite including fuel (WT in kilograms), initial mass of 
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the maneuvering fuel (W in kilograms), satellite atmospheric 
drag coefficient (CD), and satellite cross-sectional area (A 
in square kilometers). 

Motor specific impulse SPI can range from 220 seconds 
for monopropellants to 3000 seconds for ion drives. The 
default value is 230. 

The initial satellite mass WT depends on the satellite 
mission and design. A typical communications satellite has 
a mass of about 3500 kilograms. An SDI satellite may be 
considerably more massive. The default value for the 
satellite mass is 9100 kilograms. 

Initial maneuvering fuel mass W of the satellite depends 
on the design life of the satellite as well as it's 
ballistic coefficient. Low altitude satellites need at 
least one percent of their initial mass as maneuvering fuel 
if mission lives of years are to be achieved. The default 
value for initial maneuvering fuel is 100 kilograms. 

The drag coefficient CD is a measure of the slowing 
power of the atmosphere on a given satellite. In the 
thermosphere it is very pemeeenese senposiitwon dependent. A 
value of 2.2 has been settled on by orbital analysts for 
altitudes of 200 to 400 kilometers. A value of at least 3.0 
must be assigned at altitudes above 500 kilometers 
(Corresponding to the helium atmosphere). The default value 


for the drag coefficient is 2.2. 
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Satellite cross-sectional area iS again very design 


dependent. The echo series satellites had cross-sections 
measuring in hundred's of square meters. Some SDI 
satellites may also have very large cross-sections. The 


default for cross-sectional area is 50.0*107® square 
kilometers, corresponding to 50 square meters. 

The five atmosphere specific parameters are: atmosphere- 
to-earth relative rotation rate D, 10.7 cm solar flux F107 
described in Chapter II, ninety-day average solar flux FBAR, 
geomagnetic activity index AKP, and the time of vernal 
equinox passage TVE. 

The atmosphere-to-earth relative rotation rate D is the 
factor explained at the end of Chapter II. The default 
value is 0.9. 

The decimetric solar flux and the average decimetric 
solar flux F107 and FBAR are again explained in Chapter II. 
Default values of 100 are assigned to each. 

The geomagnetic index AKP is discussed in Chapter II. 
It ranges from 0 to 7 and is set at a default value of 2.00. 

With the current atmospheric model Jacchia J60, both 
FBAR and AKP are read as input but not used for the density 
computation. As seen from a sensitivity analysis, these 
parameters currently have no effect on mission life. The 
Jacchia J77 model atmosphere uses these now stranded input 
parameters. The future incorporation of the Jacchia J77 


atmosphere model will be able to use these input elements. 
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oe NAVSPASUR orbital elements are the final input 
elements in the standard input file. They are all explained 
in the orbital mechanics section, and in the code included 
in Appendix A. Briefly, they are: 

UJD: time of epoch in modified julian days 
ETU: initial mean anomaly 
GO: argument of perigee 
HO: right ascension of the ascending node 
ES: eccentricity 
XI: inclination 
AO: semi-major axis 
RND: rate of change of mean motion 
ESD: eccentricity decay rate 
ADOT: semi-major axis decay rate. 
The first seven of these are orbital elements and the final 
three are the results of perturbations. The NAVSPASUR 
orbital element set is the result of orbital measurements of 
actual satellites. 

The final block of input in the input file is the 
sensitivity analysis specification array. It must be 
included only when the sensitivity analysis mode is 
specified. The array consists of 20 lines used to make up 
the array of elements to be scanned by the program when the 
sensitivity analysis mode is selected. 

The sensitivity analysis specification array is divided 


into four columns. The first column contains the input 
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element variable name. These are the same names as 
specified in the satellite physical, atmosphere specific, 
and orbital elements input discussions. 

The second column is a control parameter. When set to 
1.0, the corresponding input element is included in the 
sensitivity analysis. Columns three and four are in a 
Fortran F16.9 format and they specify the minimum and 
maximum range of the element to be scanned. Defaults for 
the minimum and maximum have been provided in specific cases 
corresponding to logical and physical ranges of the input 
elements. In some cases zeros appear in the minimum and 
maximum columns. These zeros correspond to input elements 
not physically constrained. An example of an input element 
not physically constrained is the time of epoch (UJD). 

If a 1.0 is entered in column two of a particular row of 
the sensitivity analysis specification array, that input 
element is included in the analysis. Any number of the 
input elements may be selected for simultaneous analysis. 
When selecting a multivariable sensitivity analysis, the 
user must understand the dine e pNeReRe of the incrementing 
variables and ensure that a solution exists in the mission 
life range being considered. 

When a variable is included in the sensitivity analysis, 
the incrementing values override the initial conditions. 
All input elements not user selected for analysis (those 


with zeros in the second column) are reset to their initial 
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values after each iteration during the analysis. Appendix B 
contains sample input element files and their corresponding 
outputs. 

Input elements and output forms have been explained. 
The program flow is discussed next. The 16 subroutines and 
eight functions comprising the propellant longevity model 
are explained in Appendix A in the Fortran code. The 
routine PLEP initializes the system, including inputting 
program control parameters and initial values for all 
elements. If the sensitivity analysis mode is selected, 
subroutine SNALYS generates the array to be scanned 
linearly, complete with the necessary increments for each 
iteration. Analysis mode causes completion of all 
computations to propellant exhaustion PSIZE times. The 
routines SNALYS and PLEP are completed only once even in 
analysis mode. 

Program control is transferred from the initialization 
routine PLEP to the subroutine KOZAI, the main driver and 
bookkeeper of the system. The routine KOZAT first causes 
the NAVSPASUR Brouwer mean element set to be converted into 
a Keplerian (Kozai-like) mean orbital set, by calling the 
subroutines BRAUER, PERIOD, and MEAN. 

Subroutine MEAN calls functions XSPA-XSPM to convert the 
Brouwer osculating elements generated in BRAUER to Keplerian 


orbital elements. The Keplerian elements are then returned 
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to the subroutine KOZAI. In standard mode the mean element 
set is given as output. 

All of the above computations have simply conditioned 
the input elements to conform to the requirements of the 
propellant longevity model. The routine KOZAI now causes 
calculation of the fuel used to compensate for atmospheric 
drag after each revolution. The computation is continued 
until the fuel is exhausted or the maximum number of 
revolutions (ITER1) is reached. When the computation 
terminates, program control is transferred back to PLEP. 

To calculate the drag-compensating fuel required for 
each revolution, KOZAI calls the subroutines PERIOD, GEOP, 
THRST, and DRAG. 

The subroutine GEOP calculates the change in orbital 
elements due to the aspherical earth. It uses the J2, J3, 
and J4 terms in a Legendre polynomial expansion, as 
explained in Chapter II. Aspherical earth calculations are 
performed once each revolution. 

The subroutine THRST checks each revolution for a 
scheduled orbital adjust. If a change to the semi-major 
axis is scheduled, it calculates the change to the axis, the 
argument of perigee, the eccentricity, and the mean anomaly. 
It also computes the amount of maneuvering fuel required to 
perform the adjustment. All changes are returned to the 


subroutine KOZAI. 
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The subroutine DRAG calculates the fuel mass needed to 
overcome the atmospheric drag. This calculation is 
accomplished by calling four subroutines: SOLOC, SATLOC, 
SALT, and JAC60. 

The routine SOLOC computes the right ascension and 
declination of the sun for the model atmosphere diurnal 
density bulge. The routine SATLOC computes the right 
ascension, declination and geomagnetic latitude of the 
satellite. This is accomplished by first calling the 
subroutine POSVEL to calculate satellite position and 
velocity. The routine POSVEL in turn calls the subroutine 
NWTRPH to solve Kepler's equation (see Chapter II). The 
routine SALT computes the satellite altitude. 

Atmospheric ee is computed by the subroutine JAC60 
based on the Jacchia J60 model atmosphere discussed in 
Chapter II. A maximum, minimum, and mean density that the 
satellite will encounter in the oblate diurnal atmosphere 
are calculated. | 

The mass decrement equation (Equation 58 of Reference 1) 
is contained in the subroutine DRAG. The function MMBSIR is 
used to perform the Bessel function expansion. The routine 
DRAG returns the fuel mass decrement due to atmospheric drag 
computation to KOZAI. The routine KOZAI then decrements the 
remaining fuel by this amount and checks for fuel exhaustion 
and maximum orbital number. If fuel is less then zero, or 


the current revolution number equals ITER1, calculations are 
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terminated and control is returned to PLEP. In standard 
mode, program execution now halts. 

In sensitivity analysis mode, all input elements not 
selected for analysis are returned to their initial values, 
elements selected for analysis are linearly incremented, and 
KOZAI is again called by PLEP. The subroutine KOZAI then is 
executed, as previously described, with the new values for 
the selected elements. Upon fuel exhaustion or reaching the 
maximum revolution number, KOZAI terminates calculation and 
passes program control back to PLEP. This process continues 
until the elements user selected for analysis have been 
incremented from their minimum to maximum. When the maximum 
value of the elements user selected for analysis is reached, 
program seeeutsion terminates in the sensitivity analysis 
mode. 

The cabal output in the analysis mode is printed line- 
by-line in KOZAI upon reaching fuel exhaustion or maximum 
revolution number. The eect is a table of data, PSIZE 
long, that represents the linear scan of the specified 
elements. The Picget line is the mission life corresponding 
to the minimum value of the scanned elements. Each 
subsequent line represents the total mission life resulting 
from an increment of the elements selected for analysis and 
the initial conditions of the elements not selected. The 


final line is the mission life resulting from selected 
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elements at their maximum value and-non-selected elements at 
their initial conditions. 

Mission life is expressed in revolution number and 
modified Julian days. The last column of the tabular output 
in analysis mode represents the current value of the 
incremented input element. When more than one element is 
analyzed, the last column is the current fraction of scan, a 
real number between O and 1. 

The tabular output can easily be converted into an input 
file for a graphics routine by deleting the first 98 lines 


ef output. 


C. MODEL CODE OPERATION 

This section describes the operation of the computer 
code at the Naval Postgraduate School on the IBM 3033 
computer under VMS using the VS FORTRAN compiler. 

The propellant longevity model must be compiled in the 
VM mode using the automatic precision increase option. If 
the propellant longevity code is in filename SDI and 
filetype FORTRAN (the filetype must be FORTRAN) , the correct 
form for the compilation command is: 

FORTVS SDI (AUTODBL(DBLPAD4) ) 

The listing Senet by this compilation takes up about 20 
percent of the user's "A" disk. 

Prior to program execution, the input and output files 
must be specified. If the input file is in filename PLEPOl 


and filetype DATA2, then the correct form of the command is: 
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FILEDEF 5 DISK PLEPO1 DATA2 Al 
This informs the computer where to look for input during 
execution of the model code. If hard copy of the output is 
desired, it can be written to a file instead of the screen 
by using the FILEDEF command. If the file name is SDI and 
the desired type is DATA, the correct command is: 

FILEDEF 6 DISK SDI DATA Al 
After execution of the propellant longevity code, the 
results will be in SDI DATA on the user's "A" disk. 

To execute the model listing generated upon compilation, 
the command RUN SDI is typed (this assumes that SDI is the 
filename of the propellant longevity program). 

This program can have very long execution times. A 
general rule of thumb is that onntneien time takes less than 
PSIZE * ITER1/10,000 -Iinutes. Roughly, the computation 
takes 100 microseconds for each revolution. For example, if 
PSIZE is 100 and ITER1 is 10,000 (corresponding to a five 
year lifetime) then the IBM 3033 will take 200 mnuees to 
complete the calculations. The user must carefully choose 
small values of PSIZE and ITER1 on the initial Gnvesticatien 
of elements. 

The simplest graphical routine to construct plots of 
tabular data is EASYPLOT. The EASYPLOT routine was used to 
generate all the plots shown in Chapter IV. 

In summary, the steps of program operation are: 

(1) Compile the program using the automatic 


precision increase option. 
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(2) Select the mode of operation and input element 
values by manipulation of the input file. 


(3) Define input and output files. 
(4) Execute the program by typing the file name. 
The next chapter presents the results of the 
element sensitivity analysis performed to test the 


reasonableness during verification. 
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IV. MODEL VERIFICATION 


This chapter presents the results of the propellant 
longevity model verification at the ral Postgraduate 
School. Standard model verification tests the model against 
three questions. First, does the model address the problem? 
Second, are the model predictions reasonable? Third, how do 
the model predictions compare to real-world observations? 

This verification addresses only the first two criteria. 
Observational data for comparison with the model predictions 
is not available in an unclassified state. The quantitative 
analysis of model predictions is thus not addressed in this 
verification. 

The model does address the initial problem. It predicts 
the fuel aes decrement of a low-earth-orbit satellite doing 
intrack micro-thrusting to overcome atmospheric drag. 

The reasonableness test is the subject of the remainder 
of this chapter. An effective way to test the model 
performance is to examine model predictions resulting from 
the full “intended use" range of the input variables. 
Furthermore, if a plot is made of predicted results as a 
selected input element is incremented through its range, 
prediction trends can be found and examined. This procedure 
provides a qualitative analysis of model sensitivity to the 


selected input element. 
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To aid in these tests and trend studies, a sensitivity. 
analysis feature was incorporated into the program code of 
the model, as previously described. Graphical results of 
each single variable sensitivity analysis are presented for 
the input elements. 

The format of the following figures are all related. 
They represent plots of satellite mission life in modified 
Julian days (here after referred to simply as "days") versus 
a specified input element. All figures have the Y-axis as 
mission life and the X-axis as the variable selected for 
analysis. In each case the range of values for the input 
element is selected to remain within reasonable values of 
the input variable, as discussed in Chapter III. The 
granularity in most cases is 1/25 (PSIZE = 25). In some 
cases it is necessary to decrease the granularity to 1/100 
(PSIZE = 100) to see the actual trends. 

The default values discussed in Chapter III are used for 
all input elements not selected for analysis (except where 
noted). They correspond to a near-polar, retrograde orbit 
of an actual satellite supplied as sample input by NSWC. 
The model predicted propellant longevity using all default 
conditions and 10 kilograms of fuel is 13.74 days. When the 
default values and 100 kilograms of fuel are used, the 
resulting propellant longevity is 150.21 days. An * symbol 


on the plotted curves indicates the default value of the 
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analyzed input element. The default value names are again 
summarized here. The units can be found in Chapter III. 
Satellite physical parameters: 
SPI (motor specific impulse) 230 


WT (satellite initial mass) 9100 


W (initial fuel mass) 100 and 10 
CD (drag coefficient) 2eer and seu 
A (cross-sectional area) -00005 


Atmospheric specific parameters: 


D (atm. rotation factor) 120 

F107 (decimetric solar flux) 100 

FBAR (average F107) 100 

AKP (geomagnetic index) 220 

TVE (time of vernal equinox) 43222.7382 
NAVSPASUR orbital elements: 

UJD (time of epoch) 44619.98716775 

ETU (mean anomaly) 150.0000 


HO (R.A. of ascending node) 95.0000 


GO (argument of perigee) 200.0000 
ES (eccentricity) -0100000 
XI (inclination) 95.0000 


RND (time rate of mean motion) .000783912 
ESD (eccentricity decay rate) -.239590000 
AO (semi-major axis) 1.06 


ADOT (semi-major axis decay) -701.8652 
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The input element sensitivity analyses are divided into 
related sets. All the sets are discussed in the same basic 
format’. Each analysis is single-variable and arranged as 
follows: 
(1) The expected trend of real-world mission life 
resulting from the selected variable 
incrementing through its range. 
(2) The expected mission life trend predicted by 
the model as a result of the selected 
variable incrementing through its range. 
(3) The actual mission life results of the 
selected variable incrementing through its 


range. 


(4) A discussion of the selected variable 
sensitivity analysis. 


The explanations for the plotted trends are based on 
expected behavior in agreement with current theory. An 
understanding of the actual mechanisms producing the 
observed trends will be improved with further research. 

For a given set of satellite-orbital configurations, the 
satellite mission life will. be shortest for satellites 
traveling through the most dense atmosphere. Thus for a 
given input variable sensitivity analysis, the range of 
values corresponding to satellite travel through the most 
dense atmospheric conditions will produce the _ shortest 
mission life. The range of values of the same input 
variable corresponding to satellite travel through the least 
dense atmosphere will likewise produce the longest mission 


life. The orbital-perigee atmospheric-density-bulge 


gis 


encounter geometry is thus critical in mission life 
analysis. 

The initial sensitivity plot of an input variable is 
accomplished with 10 kilograms of fuel. In many cases, true 
geopotential-induced orbital-drift trends are not seen with 
small amounts of drag-compensating fuel. To verify trends 
caused by geopotential-induced orbital drift, the initial 
compensating fuel is increased to 100 kilograms. The 
differences in the 10 kilogram and 100 kilogram mission life 
plot trends can then be attributed to geopotential effects 
and the movement of the center of most direct solar heating 
during the mission. 

The 133.14 day nominal mission life resulting from the 
default values of the input elements and 100 kilograms of 
drag compensating fuel is almost a half-year. During this 
half-year, the earth sigs revolving around the sun causing the 
noon-sun atmospheric position to increase in longitude by 
about one degree per day. This greatly complicates trend 
analysis when using the 100 kilogram fuel case alone. [It is 
thus necessary to include the 10 kilogram fuel case. The 
13.74 nominal mission life of the 10 kilogram fuel case is 
short enough so that it shows little longitudinal movement 
of the most-direct-sun atmospheric position. 

In many of the sensitivity analysis plots, the predicted 
mission life curve has a "knee." This "knee" corresponds 


to the range of the selected input variable that produces 
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changing behavior in the mission life curve. The location 
of these ranges has significance to design research. 

The input variables are divided into four categories 
based upon their effect on mission life. The first category 
contains the input variables that affect the atmospheric- 
bulge orbital-perigee encounter geometry. The second 
category contains the input variables that affect the 
Slowing effectiveness of the atmospheric density bulge on 
the particular satellite configuration. The third category 
contains the input variables that affect the depth of 
penetration of the satellite into the atmosphere at perigee. 
Finally, the fourth category contains the input variables 
that have no affect at all on mission life. As will be 
seen, predicted mission life is most sensitive to the third 
category (factors affecting atmospheric penetration at 
perigee). | 

Factors affecting the atmospheric-density-bulge orbital- 
perigee encounter geometry are: the time of vernal equinox 
passage TVE, the time of epoch UJD, the argument of perigee 
GO, the right ascension of the ascending node HO, and the 
inclination XI. The time of vernal equinox provides a real 
time reference for earth in its orbit relative to the sun. 
When combined with the time of epoch (Sometimes called the 
time of perigee passage), these two factors determine the 
position of the noon sun in right ascension and declination 


relative to the atmosphere. The default values used in the 
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model verification position the initial noon sun at 
approximately 300 degrees longitude and 20 degrees south 
latitude relative to the celestial meridian containing the 
first point of aries. The right ascension of the ascending 
node and the argument of perigee combine to position the 
point of perigee in longitude relative to the first point of 
aries. The orbital inclination positions the point of 
perigee in latitude (or declination). Actual positioning of 
the satellite in the above orbit is accomplished by 
specifying the mean anomaly. All these orbital mechanics 
terms have been explained in Chapter II.A. 

The expected real-world behavior of satellite mission 
life as a result of varying the time of vernal equinox or 
time of epoch is discussed next. The time of vernal equinox 
TVE serves to set a real-time reference for positioning of 
the earth in its orbit relative to the sun. When the time 
of epoch UJD is specified, the position of the point of the 
most direct atmospheric heating (noon sun) in ene atmosphere 
is set. The atmospheric density-bulge follows the location 
of the most direct solar heating. The result of TVE and UJD 
varying in relation to each other is a seasonal positioning 
of the earth in its orbit around the sun at the time of 
initial perigee passage. Due to the earth's rotational-axis 
tilt, the point of most direct solar heating moves in 


longitude and latitude as the earth moves in its orbit 
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around the sun. Figures 3 and 4 in Chapter II.A aid this 
visualization. 

The earth and solar cycles effecting the atmospheric 
density~bulge combine to define its size and shape. These 
cycles are discussed in Chapter II.B, and stated again here. 
They are: diurnal, semi-annual, seasonal-latitudinal, solar 
rotational, and solar activity related. The Jacchia J60 
model atmosphere, as used in the propellant longevity model, 
incorporates only the diurnal cycle. 

Once the position, size, and shape of the atmosphere is 
set, the mission life of a drag compensating satellite 
orbiting within it can be developed. For such a satellite 
in a specific initial orbit, the following real-world trend 
in mission life (propellant longevity) should be observed. 
The trend should be cyclic, resulting from the superposition 
of all the atmospheric heating cycles. Mission lives will 
be shortest for orbits having their perigee in the most 
dense atmospheric bulges. The longer mission lives will be 
observed for satellites having their orbits farthest away 
from the density-bulge. Thus the plots of satellite mission 
life will show peaks and valleys consistent with the above 
discussion. 

In the propellant longevity model, the diurnal cycle is 
the only atmospheric heating cycle that is incorporated. 
Positioning of the earth in its orbit is accomplished by the 


combination of the time of vernal equinox and time of epoch 
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in the subroutine SOLOC. As discussed in Chapter III, the 
subroutine SOLOC positions the sun by right ascension and 
declination over the earth's atmosphere. The routine DRAG 
considers the diurnal density bulge to be centered under the 
solar position. The density profile of the atmospheric 
density-bulge is generated in the routine JAC60 based on the 
value of the input element F107. 

The propellant longevity model considers the atmosphere 
to be oblate, diurnal, and revolving around the sun with a 
23.5 degree rotational-axis tilt. The expected model- 
predicted-mission-life trend as is cyclic as TVE and UJD are 
varied. Every year should be an exact replica of every 
other year. Two peaks and two valleys in the yearly 
mission-life plot should be observed. The deepest valley 
occurs when the point of perigee is most closely centered in 
the atmospheric density bulge. With a low eccentricity (e = 
-0O1) orbit considered, a much shallower valley will occur 
about six months later. This shallower valley corresponds 
to the point of apogee being most ._closely centered in the 
atmospheric density-bulge. | 

Two peaks in plotted mission life should occur 
corresponding to the orbital pane being most perpendicular 
to the earth-sun line. The peaks should be six months 
apart, of approximately equal magnitude, and interspaced 
between the mission-life alle As the eccentricity of 


the orbit is increased, the apogee-associated valley will 
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become shallower while the perigee-associated valley becomes 
deeper. 

Figures 12, 13, 14, and 15 present the plots of the 
sensitivity analysis of TVE. In Figure 12, TVE has a range 
of 1000 days, from 38422 to 39422 days. The initial fuel is 
nO ke. The resultant predicted mission life ranges 
cyclically between 11 and 18 days. 

Figure 13 has an incremented range for TVE of 38422 to 
38522, or 100 days. Thea ec imetase lO) Kee athe resulting 
predicted mission life again ranges cyclically between 11 
and 18 days. 

In Figure 14 the width of the incremented range is 
raised to 10,000 days (about 27 years). The minimum value 
of TVE is 38422 and the maximum is 49422. Again 10 kg of 
fuel is used. The results are cyclic between 11 and 18 days 
until TVE becomes greater than time of epoch UJD. The 
results at the high end of TVE in this plot are thus 
m@eedsonabic. The time of vernal equinox passage chosen 
must occur before the satellite is launched requiring TVE to 
be less then UJD. 

The TVE is again incremented through a small range of 20 
days in Figure 15. [Initial fuel is 100 kg. The resulting 
range of mission life predicted is from 150 to 153 days. 
Note the absence of the cyclic trend. 

Figures 16 and 17 present the results of the sensitivity 


analysis plots of the time of epoch UJD with a range of 1000 
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Figure 12. Mission Life versus Time of Vernal 
Equinox (1000 Days) 
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Figure 3. Mission Life versus Time of 
Vernal Equinox (100 Days) 
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Figure 14. Mission Life versus Time of 
Vernal Equinox (10,000 Days) 
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Figure 15. Mission Life versus Time of 
Vernal Equinox (20 Days) 
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Figure 16. Mission Life versus Time of Epoch 
(1000 Days and 10 kg of Fuel) 
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Figure 17. Mission Life versus Time of Epoch 
| (100 Days and 100 kg of Fuel) 


83 


days. Note the similarity to the TVE plots, as 
theoretically expected. Figure 16 shows the 10 kg fuel 
case, and Figure 17 the 100 kg fuel case. 

The trends observed in Figures 12 through 17 "are (in 
agreement with expected model predictions. The plotted 
trends may have significant differences from real-world 
trends due to the limitations of the Jacchia J60 model 
atmosphere. It should be noted that the default values of 
TVE and UJD produce a solar location approximately 30 
degrees longitude past the winter solstice in the southern 
hemisphere heading toward the spring (vernal) equinox. The 
approximate latitude of the default solar position is 20 
degrees south. 

The right ascension of the ascending node HO positions 
the orbital plane with reference to the celestial meridian 
containing the first point of aries. When the argument of 
perigee GO and the inclination XI are specified, the point 
of perigee is positioned in longitude and latitude. This 
geometry is shown on Figure 4 in Chapter II.A. 

The real-world trend in mission life of a drag 
compensating satellite, as a result of varying HO, GO, and 
XI is cyclic. The cycle repeats every 360 degrees for HO 
and GO, and every 180 degrees for XI. The observed pattern 
of peaks and valleys will be highly variable, depending on 


the values of these three factors. 
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Model predicted mission life trends due to varying HO, 
GO, and XI through their ranges should closely parallel 
real-world expectations. This close parallel between model 
predictions and observation assumes the atmospheric density- 
bulge has been properly positioned and described. 

The default values of 95 and 200 degrees for the 
inclination and the argument of perigee combine to produce 
a point of perigee that is very nearly 20 degrees below the 
descending node in the southern hemisphere of the earth's 
atmosphere. The results of the sensitivity analysis plots 
for the right ascension of the Sacend See node and the 
argument of perigee are shown in Figures 18 through 21. 

The right ascension of the ascending node HO is 
incremented from 0 to 360 degrees. Figure 18 presents the 
10 kg fuel case and Figure 19 the 100 kg fuel case. 
Predicted mission life ranges are 11 to 18 days in Figure 
18, and 135 to 160 days in Figure 19. 

The argument of perigee GO is incremented from 0 to 360 
degrees with a 10 kg fuel case in Figure 20, and a 100 kg 
fuel case in Figure 21. The predicted mission life are 13 
to 18 days, and 148 to 160 days, respectively. 

In Figure 18 two approximately equal mission life peaks 
and valleys are observed. The valleys correspond to the 
points of perigee and apogee closely approaching the 
atmospheric density bulge. The peaks occur when the orbital 


plane is most nearly perpendicular to the earth-sun line 


85 


MISSION LIFE (DAYS) 


18 


3 14 15 


| 


i1 


mes Sak MRS i ie MN Ui i a 5) ll 





10 


0 60 120 180 240 300 360 
RA OF ASCENDING NODE (DEG.) 


Figure 18. Mission Life versus Right Ascension 
of the Ascending Node (10 kg of Fuel) 
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Figure 19. Mission Life versus Right Ascension 
of the Ascending Node (100 kg of Fuel) 
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Figure 20. 


Mission Life versus Argument of 
Perigee (10 kg of Fuel) 
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Figure 21. Mission Life versus Argument 
of Perigee (100 kg of Fuel) 
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which contains the atmospherie densicy) bullgewecen-cr The 
change in shape of the 10 and 100 kilogram curves (Figures 
18 and 19) is due to the movement of the atmospheric density 
bulge during mission life and to geopotential-induced 
OEbdtalecdc1tt. 

With default values of 95 degrees for both the right 
ascension of the ascending node and inclination, the of race 
of varying the argument of perigee from 0 to 360 degrees is 
to move the point of perigee around the initially fixed 
orbital plane (see Figure 4 in Chapter II). The plot of 
mission life corresponding to this motion should by cyclic 
showing one peak and one valley. The valley corresponds to 
the perigee point closest to the atmospheric density-bulge 
center. The peak in nilat tea mission life occurs when the 
apogee point is closest to the atmospheric density-bulge 
center. 

Figures 20 and 21 show the plotted results of the 
sensitivity analysis of the argument of perigee GO. The 
expected single peak and single valley are seen. The 
difference in the 10 kilogram case in Figure 20 and the 100 
kilogram case in Figure 21 is due to the seadonailenioueneme 
of the atmospheric density-bulge and geopotential-induced 
orbital drift. The magnitude of the mission life variations 
is less than that produced by varying the right ascension of 


the ascending node.-. This effect is due to low eccentricity 
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orbit (e =0.01) and the position of the orbital plane with 
reference to the earth-sun line. 

Factors which primarily move the orbital plane in 
relation to the atmospheric density-bulge center will show 
two peaks and two valleys in their single-cycle mission-life 
range. Factors which primarily move the point of perigee in 
an initially specified orbital plane will show one peak and 
one valley in their single-cycle range. These 
considerations assume that only low eccentricity orbits are 
being considered. 

Based on the above discussion, Figures 18 through 21 
appear to agree with expected real-world and wmodel- 
predicted mission life trends. 

The analysis of predicted mission life trends due to 
varying the inclination XI is more complex. Several factors 
must be considered in addition to the orbital-perigee 
atmospheric density-bulge encounter geometry. First, the 
diurnal atmosphere model used by the propellant longevity 
model is also oblate, and is thus more dense at the equator 
than at the poles. This fact means that satellites in near 
polar orbits travel through an overall less dense atmosphere 
than those in. a near equatorial orbit. Second, the 
atmosphere rotates in the same direction as the earth. 
Hence prograde orbital satellites experience less satellite- 
to-atmosphere . relative velocity than satellites. in 


retrograde orbits. If the satellite altitude is low enough, 
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a near-equatorial prograde orbit satellite has a longer 
mission life than the same satellite in a near-equatorial 
retrograde orbit. These additional effects complicate the 
analysis of mission life resulting from varying the orbital 
inclination. 

The default value position of the center of the 
atmospheric density-bulge is 300 degrees longitude and 20 
degrees south latitude. The bulge is moving at a rate of 
about one degree longitude per day toward the vernal 
equinox. The point of perigee is default positioned at 
approximately 275 degrees longitude and 20 degrees south 
latitude. As can be seen, the point of perigee is quite 
close to the center of the atmospheric density-bulge. The 
motion of the point of perigee, as the IMeMENSS is 
increased from near 0 degrees to near 180 degrees, describes 
a semicircle 20 anes in radius around the descending 
node. The point of perigee remains quite close to the 
initial position of the center of the atmospheric density- 
bulge. 

Figures 22, 23, 24, and 25 are the predictions of the 
sensitivity analysis plots of orbital inclination XI. The 
incremented range is from 22.5 to 157.5 degrees. Nearer 
equatorial inclinations are prohibited by the available 
arithmetic precision on the IBM 3033. Figure 22 shows the 
10 kg fuel case, Figure 23 the 100 kg fuel case, and Figure 


24 the 100 kg fuel case with a granularity of 1/100 and a 
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Figure 22. Mission Life versus Orbital 
Inclination (10 kg of Fuel) 
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Figure 23. Mission Life versus Orbital 
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Figure 24. Mission Life versus Orbital Inclination 
(100 kg of Fuel and PSIZE = 100) 
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Figure 25. Mission Life versus Orbital Inclination 
(100 kg of Fuel, PSIZE = 100, and D = 0.9) 
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Hetaewon Lactor of 0.9. Predicted mission life ranges are 
7.8 to 13.7 days in Figure 22, and 90 to 138 days in Figures 
23 and 25. 

Figure 25 shows the effect of the using 100 kilograms of 
compensating fuel and a rotation factor of 0.9. This 
situation is compared to 100 kilograms of compensating fuel 
ame@ea rotation factor of 1.0 in Figure 23. Only very slight 
differences occur between Figures 23 and 25. Mission life 
is a little longer in the retrograde orbits using a rotation 
factor of 0.9 and a little shorter in the prograde orbits 
using a rotation factor of 0.9. 

When the satellite orbit is lowered to an altitude at 
perigee of 200 kilometers (AO = 1.03) with 1000 kg of fuel, 
mission life is reduced in the retrograde-equatorial orbit 
below that of the prograde-equatorial orbit. The observed 
peaks at 70 and 140 degrees also shift leftward. The 
Seeerved dimple “seen in Foure 22 at the critical 
inclination of 63 degrees may be geopotential-induced or an 
artifact of the Kozai approximations as discussed by Taff 
[Ref. 2:pp 332-342]. 

A further explanation of the plot of mission life versus 
inclination shown in Figure 22 is presented here. A a 
detailed study of this input variable is .worthy of a 
complete paper itself. 

The low mission life found at low values of inclination 


1s due to the combination of the point of perigee position 
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and oblate-atmosphere effects. The low inclination orbits 
travel through the "fat part" of the oblate atmosphere near 
the equator. As the inclination increases, the point of 
perigee moves away from the center of the atmospheric 
density-bulge and the satellite's orbit moves up out of the 
thick equatorial atmosphere. These two factors combine to 
produce the observed rapid rise in mission life seen in 
Figure 22. Beyond near-polar values of the inclination, the 
predicted mission life again falls. The mission life valley 
seen at 120 degrees is probably due to maximum orbital time 
near the atmospheric density-bulge center. 

Figure 23 shows the effect of starting with 100 
kilograms of initial drag-compensating fuel. Note the many 
peaks and valleys. To confirm that these peaks and valleys 
represent actual predicted mission life trends and not 
Simply data scatter, the granularity was decreased to 1/100. 
Figure 24 is a plot of these results, and confirms the 
trends outlined in Figure 23. The many peaks and valleys 
superimposed over the general trend of mission life seen in 
Figure 22 are probably due to geopotential-induced orbital 
Arite. 

The results of the sensitivity and trend eerie of the 
orbital inclination XI have been discussed with reference to 
in Figures 22 through 25. Future investigation of the 


plotted data may reveal more insight toward understanding 
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the underlying mechanisms. The sensitivity analysis plots 
show no trends that appear unreasonable. 

The sensitivity analysis for the input’ variables 
effecting the orbital-perigee atmospheric density-bulge 
encounter geometry (category one) have now been presented. 
The next set of analyses are those of the second category, 
which contains the input variables affecting the slowing 
power of the atmosphere. These variables are: decimetric 
solar flux F107, drag coefficient CD, cross-sectional area 
A, specific impulse SPI, and the atmospheric rotation factor 
D. The initial mass of drag-compensating fuel is also 
included in this category because it improves satellite 
mission-life resistance to atmospheric drag. 

The real-world expected effect on mission life as a 
result of varying the decimetric solar flux F107 is to vary 
the heating effectiveness of the sun. As the decimetric 
solar flux increases, the density of the atmosphere will 
increase. 

The Jacchia J60 model atmosphere used in the propellant 
longevity model subroutine JAC60 uses the F107 flux as a 
multiplier for the exponential density term. The expected 
model predictions of mission life as the F107 flux is varied 
through its range of 50 to 300 is a multiplicatively 
decreasing plot. 

The results of the sensitivity analysis for predicted 


mission life as the decimetric solar flux is varied are 
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presented in Figures 26 and 27. Figure 26 is the result of 
10 kg of initial maneuvering fuel, and Figure 27 is the 
result with 100 kg of fuel. The predicted mission life in 
Figure 26 ranges from 24.99 to 3.96 days. The 100 kg fuel 
case in Figure 27 predicts a mission life ranging from 
270.12 to 43.07 days. These results are in agreement with 
expectations. 

The next set of analyses consider effects of the drag 
coefficient, satellite cross-sectional area, and specific 
impulse. All of these factors are multipliers acting to 
change the resistance of the satellite to atmospheric drag. 
The drag coefficient CD is a proportionality constant in the 
drag force model. The satellite cross-sectional area A 
changes the size of the satellite that the atmosphere 
actually acts against. As satellite cross-sectional area 
increases it will take multiplicatively more energy to 
compensate for the orbital energy lost by the satellite 
working against the drag force. The specific impulse SPI is 
a multiplier that expresses the effectiveness of the 
expelled fuel mass in overcoming the drag force. 

The rotation factor D is the ratio of the atmospheric 
rotation rate to the earth's rotation rate. Its effect is 
to change the satellite to atmospheric relative velocity as 
the inclination changes. A logical range of D is between 
the values 0.9 and 1.0. The sensitivity analysis range of 


0.1 to 10.0 was chosen to show model sensitivity for a wide 
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range of D. As presented in Chapter II Equation 12, the D 
factor is combined with the inclination and the ratio of the 
radius-at-perigee to the velocity-at-perigee producing the 
§ factor. This factor is used to convert the satellite 
tangential velocity into the  satellite-to-atmospheric 
relative velocity. 

The propellant longevity model sen AL CD, 9S, and SPI 
aS proportionality constants, multiplicatively changing the 
effectiveness of the expelled fuel mass. The routine PLEP 
combines these four factors into a single constant DEL given 


by: 
LEE = A> Cb = Sues. 14159/SPL (14) 


This constant DEL is passed to the routine KOZAI which in 
turn passes it to the routine DRAG. In the routine DRAG, 
DEL is used as the final multiplier in determining required 
fuel for each orbit. 

The expected model-predicted mission-life as a result of 
changing these factors is a multiplicative increase when A, 
CD, or 6 are increased linearly. A multiplicative decrease 
should occur when SPI is increased linearly. 

Expected results of varying D from 0.1 to 10.0 in the 
near-polar retrograde orbit are shown in Figure 28. Note 
that a very slight lowering of mission life is seen when the 


rotation factor increases. Figures 29 and 30 show the 
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plotted result of varying D through the same range of 0.1 to 
10.0. However, in Figure 29 an inclination of 25 degrees is 
used instead of the default 95 degrees, and in Figure 30 an 
inclination of 155 degrees is used. In Figure 29 it is seen 
that the mission life multiplicatively increases as D is 
linearly increased. Figure 30 shows a multiplicative 
decrease in mission life as D is linearly increased. These 
trends are due to the change in the satellite-to-atmosphere 
relative velocity. The increasing mission life expectancy 
within a strongly prograde orbit (Figure 29) is a result of 
the decreasing satellite-to-atmosphere relative velocity. 
The decreasing mission life expectancy within a strongly 
retrograde orbit (Figure 30) is a result of increasing 
satellite-to-atmosphere relative velocity. The plots shown 
in Figures 28, 29 and 30 are in agreement with expected 
results. 

The results of the sensitivity analysis plots on mission 
life as a result of varying CD, A, and SPI are shown in 
Figures 31 through 35. 

Figures 31 and 32 result from varying the satellite 
drag coefficient CD from 1.0 to 3.5. Figure 31 represents 
the 10 kg fuel case and Figure 32 the 100 kg case. The 
predicted mission life range in Figure 31 is 31.18 to 8.63 
days, and in Figure 32 from 308 to 99.58 days. 

The sensitivity analysis mission-life trend plots of 


satellite cross-sectional area A are shown in Figures 33 and 
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Coefficient (100 kg of Fuel) 


AE oh, 


MISSION LIFE (DAYS) 


80 


60 


00 


40 


30 


Se 10 KG OF FUEL 





0.0 1.0 | 2.0 3.0 
CROSS-SECTION (SQ. KILOMETERS) *10- 


Figure 33. Mission Life versus Cross-Sectional 
Area (10 kg of Fuel) 
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Figure 34. Mission Life versus Cross-Sectional 
Area (100 kg of Fuel) 
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Figure 35. Mission Life versus Specific 
Impulse (10 kg of Fuel) 
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34. Figure 33 is the 10 kg fuel case and Figure 34 the 100 
kg case. In both cases, A was varied from 0.000010 to 
0.0001000 square kilometers. The ranges mapped by these 
plots are 80 to 0.64 days and 736.29 to 108.7 days, 
respectively. 

Figure 35 is the analysis of motor specific impulse SPI 
varying from 200 to 3000 with 10 kg of initial fuel. The 
apparent linear shape of the SPI plot is due to its inverse 
multiplicative action and the range of the constant DEL (14) 
over which it acts. As SPI increases, the drag compensating 
fuel is more effective. Linearly increasing SPI has the 
same effect as starting at the default value * on the cross- 
sectional plots and moving along the plot to the left. Note 
oat the shape of cross-sectional mission life plot in this 
range is closer to the expected mirror image of the SPI 
DiLGiES 

As noted earlier, mission life plots showing "knees" 
have eWentieicance. to research and design studies. The knees 
in the cross-sectional area mission life plots are very 
pronounced. Less distinctive knees are seen in the mission- 
life plots of the decimetric solar flux and the drag 
coefficient. The default value of the cross-sectional area 
is right in the knee range of the predicted mission-life 
plots, as seen in Figures 34 and 35. 

The final input variable in the second category is the 


initial propellant mass W. The expected mission life of a 
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drag compensating satellite as fuel is increased linearly 
should be a linear increase in the predicted mission life. 
A given amount of fuel produces a certain mission life. 
Twice that amount of fuel in the same general configuration 
should produce about twice the mission life. The expected 
shape of the sensitivity analysis plot of predicted mission 
life for linearly increasing initial propellant mass is a 
approximately straight line. Figure 36 is the result of 
varying W on mission life. The plot is approximates a 
straight line as expected. The ripples are probably due to 
geopotential-induced orbital-drift or to movement of the 
atmospheric density-bulge during mission life. 

The results of the sensitivity and trend analysis 
mission-life plots for CD, A, D, SPI, and Ww agree with 
expectations. All input variables in the second category 
have been examined. The next input variables to be examined 
are in the third category. 

Consider now the input variables that affecting the 
depth of penetration of the satellite into the atmosphere at 
perigee. The variables that affect the satellite altitude. 
at perigee are the semi-major axis and the orbital 
eccentricity. These variables show the most pronounced 
effects on mission life. The real-world effect on mission 
life when linearly changing the altitude at perigee is an 


exponential increase in mission life. The effect of varying 
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Initial Propellant Mass 
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the eccentricity and the semi-major axis on the altitude at 


perigee is summarized by 
Ipn=a* Gite e) (15) 


Here, Yp is the radius at perigee, a the semi-major axis, 
and e the eccentricity. A linear increase in the semi- 
major axis results in a modified exponential increase in 
mission life. A linear increase of the eccentricity results 
in a modified exponential decrease in mission life. 

The propellant longevity model computes the satellite 
altitude at perigee in the subroutine SALT. The subroutine 
DRAG then passes the altitude at perigee to the routine 
JAC60 where it is used to develop an exponent for an 
empirical model of atmospheric density. 

Mission life results of varying the orbital eccentricity 
ES from 0.01 to 0.1 are presented in Figures 37 and 38. 
Figure 37 is the 10 kg fuel case, and Figure 38 the 100 kg 
fuel case. The predicted mission life ranges from 14 to 0 
days in Figure 37, and 150 to 0O days in Figure 38. 

| An analysis of the semi-major axis AO is presented in 
Figure 39 with a range of 1.025 to 1.175 earth radii. In 
this case, 1 kg of fuel is allocated initially. The 
resulting mission life is from 0 to 384 days. 

The plots of mission life resulting from varying the 


semi-major axis AO and the orbital eccentricity ES, show 
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Figure 38. Mission Life versus Orbit 
Eccentricity (100 kg of Fuel) 
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Figure 39. Mission Life versus Semi-Major 
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major knees in their curves. The * symbol indicates the 
default value of the varied element. Note that the amount 
of fuel used in the analysis of the semi-major axis is 1 kg. 
The use of the usual 10 kg produces a plot that exceeding 
the reasonable range of mission life. The results shown in 
Figures 37 through 39 are in agreement with expectations. 

The input variables in the fourth category (no effect on 
mission life) are the mean anomaly ETU, initial satellite 
weight WT, average decimetric solar flux FBAR, and the 
geomagnetic index AKP. FBAR and AKP have no effect on 
mission life due to the limitations of the Jacchia J60 
atmosphere model. In the real world these two factors may 
combine with the decimetric solar flux F107 to puoatee 
changes in the atmospheric density of three orders of 
magnitude. 

The model theory used by Dr. Parks [Ref. 1] to predict 
the propellant longevity of a satellite doing intrack micro- 
thrusting to overcome atmospheric drag is independent of the 
satellite mass. The fuel required to overcome drag-induced 
orbital energy loss is based strictly on the satellite's 
shape and size. 

Figure 40 presents the sensitivity analysis plot of 
predicted mission life as a result of varying the satellite 
initial mass WT. The plot is a straight horizontal line, as 


expected. 
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The final analysis is concerned with the mean anomaly. 
In the near-circular default orbit, the mean anomaly and the 
eccentric anomaly are nearly the same. The Keplerian mean 
anomaly functions to position the satellite in its 
designated orbit. For a mission life lasting several days, 
it should make no difference where the satellite is 
initially positioned in its orbit. Rather, it's the orbit 
itself that makes the difference. 

Analysis plots of the mean anomaly ETU are presented in 
Figures 41 and 42. The incremented range is from 0 to 179 
degrees. Figure 41 is the 10 kg fuel case and Figure 42 is 
the 100 kg fuel case. Predicted results range from 13.74 to 
13.87 days in Figure 41, and from 149 to 154 days in Figure 
42. 

The results in Figure 41 are as expected, an approximate 
horizontal line. In Figure 42, a valley is seen in the 40 
to 120 degree range. The valley is only four days deep, but 
is difficult to explain. Further research may reveal its 
underlying mechanism. | 

The results of the trend and sensitivity analysis have 
been presented. Mission life has been found to be most 
sensitive to those factors that affecting the orbital radius 
at perigee. The main unexpected results are found in the 
analysis of the inclination and the mean anomaly. The 
pattern of many peaks and valleys found in the 100 kilogram 


fuel case (Figures 23, 24 and 25) warrants further research. 
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Figure 41. Mission Life versus Mean 
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Anomaly (100 kg of Fuel) 
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The four day dimple seen in the mean anomaly plot (Figure 


41) also warrants further study. 
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V. SUMMARY 


This final chapter presents the summary and conclusions 
of the thesis. Recommendations for model use are also 
given. 

The low-earth-orbit satellite propellant longevity model 
has been tested in a sensitivity analysis. Based on the 
results of the sensitivity analysis the model reasonableness 
is confirmed. The first chapter presented a model overview, 
model motivation, and the process employed to verify the 
model. Chapter I also discusses some changes made to the 
model computer program for the IBM 3033 at the Naval 
Postgraduate School. Chapter II presents the orbital 
mechanics a8 atmospheric concepts as background for model 
understanding and it's operation. In Chapter III the model 
and its computer program were discussed. Operating 
instructions specific for the computer program use at the 
Naval Postgraduate School are also presented there. The 
sensitivity analysis results are presented in Chapter IV. 

The propellant longevity. model predicts how long the 
fuel of a low-earth-orbital satellite thrusting to overcome 
atmospheric drag will last before starting uncontrolled 
decay into the atmosphere. Model verification reasonable- 
ness tests reveal no inconsistencies with the established 


theory. The results of the sensitivity analysis are 
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therefore in Be eemnent with the expectations. The 
irregularities in the model predictions (occurring with 
inclinations of less than 22.5 degrees) are due to computer 
round off error, not the model itself. 

The propellant longevity model currently uses’ the 
Jacchia J60 model atmosphere. This atmosphere model may 
Significantly reduce the real-world accuracy of predicted 
mission lifetimes. 

The short term confidence factor in predicted mission 
life might be increased by the employment of the Jacchia J77 
model atmosphere. One drawback to the use of the Jacchia 
J77 model atmosphere is that it may double the computer run 
time. 

One solution to the trade-off of increased accuracy of 
the Jacchia J77 model versus the shorter esineeevehon time of 
the Jacchia J60 model, would be to allow the user to select 
an atmosphere option. 

Transportation of the computer program of the model to 
an IBM AT could provide more flexibility in preliminary 
design analysis. 

This model is a tool for the low-earth-orbit satellite 
system planner in deSigning and managing satellite systems 
that have precise orbital element control requirements. 
These systems are those with SDI and reconnaissance 


missions. 
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This model, with the addition of the Jacchia J77 model 
atmosphere and a subroutine to calculate fuel required to 
compensate for aspherical geopotential forces, could be used 
to optimize SDI defensive satellite constellation design. 

Existing real-world data to compare with the model 
predictions is limited to a very classified system. The 
final step in verifying the propellant longevity model would 


be to make this comparison test. 
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PSE CARO Se ank So, Nasco UGeWAs DEEN TRANorORTED-1T0 THE IBM 3033 
AND TRANSLATED TO FORTRAN 77 BY LT. C.D. NOBLE. ADDITIONAL DOC- 
UMENTATION WAS ADDED DURING THE TRANSPORTATION AND TRANSLATION. 
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ARE SITED BY NUMBER. THEY ARE: 


REFERENCE 1 BROUWER, D., “SOLUTION OF THE PROBLEM OF 
ARTIFICIAL SATELLITE THEORY WITHOUT DRAG", 
THE ASTRONOMICAL JOURNAL, VOL 69%, NO 1274, 
Big) Oa 59 tet 5D . 


REFERENCE 2 AGRAWAL, B.N., "DESIGN OF GEOSYNCHRONOUS 
SPAGEGRAFI", FRENTICE-HACE, 1986. 
REBERENCE 3 JACCHIA, L.G., "THERMOSPHERIC TEMPERATURE, 


DENSITY AND COMPOSITION: NEW MODELS", 
SPECIAL REPORT 375, SMITHSONIAN ASTRO- 
PHYSICAL OBSERVATORY, 1977. 


REFERENCE 4 JACCHIA, L.G., "EMPIRICAL MODELS OF THE 
THERMOSPHERE AND REQUIREMENTS FOR IMP- 
a adl AFGL~TR-81-0195, ADA101946, 
cog 


REFERENCE 5 TACCHI AV umG. > meANALYSIS OF DATA FOR THE 
DEVELOPMENT OF DENSITY AND COMPOSITION 
MODELS Or IHe UPPER ATMOSPHERE", AFGL-TR- 
81-0230, AD1064920, 1981. 


REFERENCE 6 ORRsebeHen VORBITAL LIFETIME PROGRAM", 
NASA TM 87587, 1985. 


REFERENCE 7 TAFF, L.G., "CELESTIAL MECHANICS: A COMP- 
UTATIONAL GUIDE FOR THE PRACTIONER", 
JOHN WILEY & SONS. INC., 1985. 


REEERENCE $ FUGAGEE, R.G., “AN INTRODUCTION TO ATMOS- 
BHWeXxLCePhyolCs", “ACADEMIC PRESS, 1930. 
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MODULE FUNCTION 
Pilee SYSTEM INITIALIZATION AND CONTROL. 
CALLED BY: CALLS: 
CSUBROUTINES) TIMES CALLED 
KNX SNALYS 1 
KOZAI 1/ ITERATION 
SNALYS SENS LLL MODE SCANNED ELEMENT GENERATOR. 
CALLED BY CALLS: 
CSUBROUTINES) TIMES CALLED 
PLEP KKK KKH 
KOZAI SYSTEM WORKHORSE, DRIVER AND RECORD KEEPER. 
NSWC TR 83-135. 
GACUED BY CALEES: 
CSUBROUTINES) TIMES CALLED 
PEER BRAUER l 
PERIOD 1 
MEAN l 
THRST 1/CSAT.REV.) 
P GEOP I7CSAT SREY 2) 
DRAG PACSAT, REVS® 
BRAUER CONVERTS NAVSPASUR BROUWER MEAN ORBITAL ELEMENT 


SET INTO A BROUWER OSCULATING ORBITAL SET. 
A CODING OF THE BROUWER-LYDDANE THEORY. 
NSWC TR 83-107, 83-155, BROUWER CREF. 1). 


CALLED BY: CALLS: 
(SUBROUTINES) TIMES CALLED 
KOZAI KK KNX 
POSVEL COMPUTES SATELLITES POSITION AND VELOCITY IN 


INERTIAL SPACE FROM THE OSCULATING ORBITAL 
ELEMENTS. NSWC TR 82=387 (PP 2G—25" 


CALLED BY: CAGES: 
(SUBROUTINES) TIMES CALLED 
SATEUG NWTRPH 1 
NWTRPH A NEWTON-RAPHSON SOQUUTION TO KEPLERS EQUATION. 
NSWC TR 82 387 P 23. 
CALCEED=5r-: CALLS: 
(SUBROUTINES) TIMES CALLED 
POSVEL KKK Pees 3 
FUNCTIONS COMPUTES SHORT PERIODIC EFFECTS ON ORBITAL 
XSPA-XSPR ELEMENTS. CALLED BY SUBROUTINE MEAN. NSWC 81-456 
PAGES 9-11 AND NSWC TR 83 243 REF 3. 
CALLED BY: CA IeS: 
CSUBROUTINES) TIMES CALLED 
MEAN XK KKK 
MEAN USES THE WALTER METHOD TO CONVERT OSCULATING 


ORBITAL ELEMENTS TO KOZAI LIKE MEAN ELEMENTS 
USING FIRST ORDER PERIODIC VARIATIONS. 
SIMILAR TO THE 4ETHOD TN NSWE TR 83 135 PP I6-18- 


CALLED BY: CALLS: 
CFUNCTIONS) TIMES CALLED 
KOZAT XSPA-XSPR 1 
GEOP COMPUTES AVERAGED SECULAR CHANGES IN KOZAI MEAN 


ELEMENTS INDUCED BY GRAVITATIONAL FIELD THROUGH 
FOURTH ZONAL HARMONIC (J2-J4 TERMS). NSWC TR 81- 
456. “EGS ce .5)'. 


CALLED BY: CALLS: 
(SUBROUTINES) TIMES CALLED 
KOZATI Pes S 3 ; Pes = 3 
SOLOC COMPUTES THE RIGHT ASCENSION AND THE DECLINATION 
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OF THE SUN FOR THE MODEL ATMOSPHERE. NSWC TR 81- 
BOO mne.  2O—ei., 


CAELED BY: CALLS: 
CSUBROUTINES) TIMES CALLED 
DRAG XXX KKK 
PERIOD COMPUTES THE ANOMALISTIC PERIOD. 
CALLED BY: CALLS: 
CSUBROUTINES) TIMES CALLED 
KOZAI KX KX % XX 
DRAG COMPUTES DRAG EFFECTS AND FUEL MASS DECREMENT 


REQUIRED TO MAINTAIN THE SEMIMAJOR AXIS VIA 
IN TRACK MICROTHRUSTING. NSWC TR 82-243 EQ. 
58. THIS ONE EQUATION IS THE HEART OF THE MODEL. 
CALLED BY: CALLS: 
CSUBROUTINES) TIMES CALLED 
KOZAI SOLOC 


JAC60 

- CFUNCTIONS) 
MMBSIR 
XSPM 


JAC60 COMPUTES ATMOSPHERIC DENSITY BASED ON THE JACCHIA 
(60) MODEL ATMOSPHERE. CALLED THREE TIMES TO 
CONCUIC OULEFERENT VALUES FOR THE DENSITY FOR EACH 
REVOLUTION: MAX, MIN AND OP. 


—r NA 


CALEED BY: CARLES: 
CSUBROUTINES) TIMES CALLED 
DRAG KX KKK 
THRST CALCULATES CHANGES TO KOZAI ORBITAL ELEMENTS IF 


AN ORBITAL ADJUST IS SCHEDULED. CALCULATES THE 
MASS OF FUEL USED TO PERFORM THE ADJUST. 
NSHC TR 83-31 EQ. 16, AND NSWC TR 81-456 EQ. 2-26 


CAPPED BY: CALLS: 
CSUBROUTINES) TIMES CALLED 
KOZATI XXX KX XX 
SALT ae a SATELLITE ALTITUDE. NSNC TR 81-456 EQ. 
GALLED BY: CALLS: 
CSUBROUTINES) TIMES CALLED 
DRAG KK KX KX 
SATLOC COMPUTES SATELLITE LOCATION. RIGHT ASCENSION, 
DECLINATION AND ITS GEOMAGNETIC LATITUDE. 
CALLED -6Y: CALLS: 
é (SUBROUTINES) TIMES CALLED 
DRAG POSVEL 1 
FUNCTION COMPUTES BESSEL FUNCTIONS. AN IMSL ROUTINE 
(MMBSIR) RESIDENT ON THE IBM 3800. THIS WILL VARY WITH 
EOGAL INSWAREATIONS. 
GALLED* BY: GALES: 
(SUBROUTINES) TIMES CALLED 
DRAG KRX¥ XXX 
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PROGRAM PLEP xX 
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geo 5 


PROGRAM PLEP 


NOA 


She 
WT 


W 
CD 


ATMOSPHERIC 
D 


THIS MODULE IS THE INPUT, 


IRVCK) 


DACK) 


PSEZE 


ars 


DESCRIPTION OF ARGUMENTS 


INPUT ARGUMENTS 
PROGRAM CONTROLS 
ITERL__ 
ITER2_____ 


INTEGER, MAXIMUM NUMBER OF REVOLUTIONS DESIRED 
FOR CONSIDERATION. A CONTROL PARAMETER. 

INTEGER, ORBITAL PRINT FREQUENCY. A CONTROL 
PARAMATER. SETTING ITER2 = TO 100 WILL REDUCE 
Uabsecte OUTPUT BY ONLY PRINTING OUTPUT EVERY 100 
ORBITS. 


INTEGER, CONTROL PARAMETER ee oense) SPECTREIES 


SENSITIVITY ANALYAIS MODE. PSENS=0 SPECIFIES 
STANDARD PROPELLENT LONGIVITY MODE. 
INTEGER, CONTROL PARAMETER. IET=1 PRINTS THE 


MHESSPECTPIC SARELLITE 


HEH HHH HH HEH HH HH HAH HH HK HH HHH HH HHH HHH HH HHH HK HHH HHH HHA HH HH HH HHH HK KKK KK KKK 
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INITIAL ANOMALISTIC PERIOD. IET=0 SUPRESSES PRINT.% 
SYSTEM RESET TO 0 WHEN SENSTRIVityY MODE 15) SEEECT ED 


INTEGER, NUMBER OF SCHEDULED ORBIT ADJUSTS. THIS 
IS A CONTROL PARAMETER WITH ALLOWED VALUES FROM 
ZERO TO TEN. THE ONLY ORBIT ADJUSTS CURRENTLY 
SUPPORTED ARE CHANGES TO THE SEMI-MAJOR AXIS. 
INTEGER ARRAY. (K=10 MAX).THE ELEMENTS OF THE 
ARRAY SPECIFY THE REVOLUTION NUMBER DURING WHICH 


EACH OF THE ORBITAL ADJUSTS ARE PERFORMED. K MUST 


EQUAL NOA. 

INTEGER ARRAY. €K=10 MAX) THE ELEMENTS OF THE 
ARRAY SPECIFY THE CHANGE IN KILOMETERS OF THE 
SEMIMAJOR AXIS. K MUST EQUAL NOA 


INTEGER, CONTROL PARAMETER THE NUMBER OF INTERVALS 


IN. THE SENSITIVITY ANAEYSIS: 


SATELLITE PHYSICAL PARAMETERS 


REAL, MOTOR SPECIFIC IMPULSE IN SECONDS. AGRAWAL 
CREF. 2:PP. 163-166) GIVES Sau2Ue VALUES: 

REAL, INITIAL MASS OF THE SATELLITE IN KILOGRAMS. 
CINCLUDING THE FUEL MASS) 

REAL, INITIAL MASS OF FUEL IN KILOGRAMS. 

REAL, DRAG COEFFICIENT OF Jade Sane eer ie. 
DIMENSTONLESS. NORMALLY VARIES FROM ONE TO 
AROUND THREE, DEPENDING ON VELOCITY, DENSITY, 
COMPOSITION OF ATHOSPRERE ANDeEar SICAL DESTON. 
IS MOST OFTEN SET AT 2220 S39R AETIRUDES SEINEER 
200 AND 400 KM. FOR ALTITUDES OF 1000 «M CD [5S 
SET AT JUST ABOVE 3.0. JACGHITA CREF Ss: FaZe 
REAL, SATELLITE CROSSESTIONAL AREA IN SQUARE 
KIMOMETERS. 


SPECIFIC PARAMETERS 
REAL, RATIO OF ATMOSPHERIC TO EARTH ANGULAR 


ROTATION RATES. A CONSTANT IN EQUATION 11 OF 
NSWC TR 83-243. UNITLESS. THE ACTUAL VALUE‘WILL 


DEPEND ON LATTITUDE AND ALTITUDE. RANGE .8 TO 1.0. 
REAL, DECIMETRIC SOLAR FLUX. JACCHIA CREF 4 & 5). 


THIS VALUE IS USED BY THE INDUSTRY TO REPRESENT 
THE XUV FLUX IN THE THERMOSPHERE. THE XUV FLUX 
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HE HEHE HE HE IE HE HE HE HE HE HE HE 3K DE HE 3 3 3 DE DE HE DE 3 DE HE HE DE 3 HE HE HE DE DE HE HE HE IE HE HE EE DE DE HE HE DE DE DE 2K DE IEE DE DE EEE DE DE DE HE DE HE EE DE OE HEE 
CONTROL AND INITIALIZATION SECTION FOR THE 
REST OF THE SYSTEM. IT ACCEPTS THE 'FIVE CARD' FORMAT OF NAVSPASUR 
AS THE BROUNER MEAN ORIBITAL ELEMENT SET. 
PARAMETERS ARE ENTERED. ATMOSPHERIC MODEL ELEMENTS ARE INITIALIZED. 
FINALLY, SYSTEM CONTROL CONSTANTS ARE SET TO GENERATE THE OUTPUT. A 
FURTHER DESCRIPTION OF THE ABOVE FUNCTIONS IS CONTAINED IN THE 
DESCRIPTION OF ARGUMENTS. 

PESTS SSS ESET SES ESET SSS SES PESTS SSE SESS SELES ST SSS ESET ESSE SST SESE S STS SS TELS S| 


«K 


% 
% 
X% 
% 
% 
x 
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FBAR 


AKP 


TVE 


NAVSPASUR ORBITAL ELEMENTS 


THE NAVSPASUR ORBITAL ELEMENT SET IS THE BROUWER MEAN 
BECMenNy sel. SOROUNCK \Ser. L) tome DES) DESCREPTIGN FOR 
PAS ELEMENTS SET. 


UJD 


ETU 
G0 
HO 


ES 


XI 


RND 


ESD 
AQ 


IS THE FACTOR CREATING THE SUN FACING DIURNAL 
DENSITY AND TEMPERATURE BULGE IN THE THERMOSPHERE. 
IT IS DIRECTLY RELATED TO SOLAR ACTIVITY. IT 
VARIES WITH THE 27 DAY SLOAR ROTATION AND 11 YEAR 
SUNSPOT ACTIVITY CYCLE. ORR CREF. 6:APPENDIX B) 
INDICATES RANGES OF 73.3 TO 242.9. XUV FLUX AND 
THE GEOMAGNETIC HEATING EFFECTS CAN CAUSE 

DENSITY TO CHANGE BY THREE ORDERS OF MAGNITUDE AT 
SOS KTROMETERS. THE V0.7 CM FLUX IS DEGRABLE FROM 
THE GROUND AND THEORETICALLY PARELLELS THE XUV 
FLUX. THE XUV FLUX IS NOT MEASUREABLE FROM THE 
GROUND DUE TO NEAR COMPLETE UPPER ATMOSPHERE 
ABSORPTION. 


mK KK MK MK OK MK MK OK OK OM OM 


REAL, THE AVERAGE 10.7 CM FLUX. THE EFFECTS OF THEX 


XUV HEATING DO NOT RETURN TO THE ORIGINAL UNHEATEDX 
CONDITION IN AN EARTH ROTATION. A HIGH HEATING 
CONDITION WILL HAVE SOME LONGER TERM EFFECTS. 
VHESe SGUGER Term EFFECTS ARE MODELED HERE BY 
TAKING A MOVING 90 DAY 10.7 CM AVERAGE. 

FBAR RANGES FROM A LOW OF 73 TO A HIGH OF 230. 
REAL, GEOMAGNETIC ACTIVITY INDEX, JACCHIA (CREF. 4% 
AND 5). GEOMAGNETIC ACTIVITY ALSO CAUSES HEATING 
IN THE THERMOSPHERE. THE EXACT MECHANISM IS NOT 
UNDERSTOOD. IDEALLY QUEIT MAGNETIC CONDITIONS 
CORRESPOND TO AN INDEX OF 0. MAXIMUM ACTIVITY FOR 
WHICH SIGNIFICANT DATA EXISTS ARE RELATED TO AN 
INDEX OF BETWEEN 6 AND 7. NORMAL ACTIVITY 
INDICATES AN INDEX OF BETWEEN 1.AND 2. 
REAL, TIME OF VERNAL EQUINOX PASSAGE. USED TO 
TRANSFORM TIME FROM INERTIAL TO EARTH FIXED 
REFERENCE FRAMES. FOUND IN THE ASTRONOMICAL ALMANA 
FOR THE YEAR OF CHOICE. (MODIFIED JULIAN DAYS). 
ASDESCRIPTIONSIS SIN TAFF (REF. 7:PP. 105-104). 
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REAL, TIME OF EPOCH IN MODIFIED JULIAN DAYS. 
THIS IS THE TO IN THE EQUATION FOR THE MEAN ANOMALY 
oeente=ett0 (tah eCh = 90) ). TAFF CREF. 7). 
REAL, INITIAL MEAN ANOMALY MO. MEASURED IN DEGREES* 
HAS A RANGE OF 0 TO 360. 

REAL, MEAN ARGUMENT OF PERIGEE. RANGES FROM 0 TO 
360 DEGREES. NOT DEFINED FOR CIRCULAR ORBITS. - 
REAL, MEAN RIGHT ASCENSION OF THE ASCENDING NODE. 
RANGES FROM 0 TO 360 DEGREES. NOT DEFINED FOR 
EQUITORIAL ORBITS 

REAL, MEAN ECCENTRIY. FOR THIS THEORY RANGES FROM 
Omer 1000.1. THESE ARE LOW ECCENTRICITY ORBITS. 


REAL, MEAN INCLINATION. RANGES FROM 0 TO 180 
DEGREES. VALUES OF 0 AND 180 ARE IN EQUITORIAL 
ORBITS. INCLINATION OF 90 IS A POLAR ORBIT. RANGES: 
OF 0 TO 90 ARE IN PROGRADE ORBITS. INCLINATIONS 
90 TO 180 ARE IN RETROGRADE ORBITS. FOR Atl 
INCLINATION OF 0 OR 180 HO IS NOT DEFINED. THERE 
IS NO ASCENDING NODE IN EQUITORIAL ORBITS. 

REAL, RATE OF CHANGE OF MEAN MOTION. MEASURED IN 
REVOLUTIONS PER CDAY SQUARED) .VALUES OF .00001 
AND ABOVE ARE CONSIDERED SIGNIFICANT WHEN 
ATMOSPHERE INDUCED. REFERENCE (TAFF) at SY / 
REAL, ECCENTRICITY DECAY RATE MEASURED IN 
DHVERSESSECONDS . 

REAL, KAULA SEMIMAJOR AXIS. MEASURED IN 

MEAN EQUITORIAL EARTH RADII (6378.165 KM). RANGES * 
FROM 1.02594 TO 1.16847 TO MAINTAIN AN ALTITUDE OFX 
100 TO 1000 KM FOR RADIUS AT PERIGEE. THIS MUST BEX 
CAREFUL L Se SCALED TO ECCENTRICITY TO PREVENT 5 
CRASHING THE SATELLITE. PERIGEE ALTITUDES ABOVE x 
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600 KM ARE ESCENTIALLY AVOVE THE ATMOSPHERE. 
ADOT REAL, SEMIMAJOR AXIS DECAY RATE. 


x 
x 
x 
x 
x 
x 
x 
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222 22 Oe ttt tee Peete tee PPS SSeS OTE SS SESS SES SEES SSS SSE LES ESET Te F | 


DIMENSION SNSARC20), SNSMINC20), IRVC10), DAC1O), SNSDELC20) 


2.22.2, 0,0,2,2, 2, 0,0,%,.0,09,2,.2,0,2,.0,0,2.0,0,.0,2,2,0,2,2,.2, 2,2, 0, 5,050,050 2 2 St et eet steseseseseses es & ; 


* THE SET OADJ IS COMMON ‘TO PLEP, KOZAI, THRST AND SNALYS. x 
HEH BE HE DE DE EE DE HEE BE HE BE HE DE DE HE HE DE DE HE DE HE DE DE DE DE DE HE HE DE DEE HE DE IE OC DE DE HEE HE IE HE HE HE DE HE HE HEH HEHE HE HE HE HE HEHEHE HE EE EK NEN 


COMMON 7 OADJ/ NOA, IRV, DA, SPI, WT, PSENS 
X 
PES eSESS EEE SEE Sse SS SeCe SCE Sree ee ete eee rere ete eee eee eee ee eres es eset eres: 


ES THE SET XXXX IS COMMON TO PLEP AND PERIOD. x 
HE DE 3 HE DE 3 HE HE DE DE HE HE DE DE HE HE DE DE HE HE DE DE HE DE DE DE IE HE DE DE HE HE HE I HE HE DE DE HE HEE DE DE HE HE HEH HE HE HE HE DE HE HE HEH HEHE HE HEHE HE HE HE HE HE 


3 
COMMON 47 XXXX / IET 
x 
HE HEHE DE HE HE HE DE HE HE HC HEHE DE HE HE DE HE DE HE HE HEHE DE HE HE HE HE HE DE IE DE HE HE HE DE HE DE DE DE HEE HE DE DE DE DEE HE HE DE HE DE DE DE HE HE HE HE DE DE EDK TE HE HE HE HE 


* THE SET SOLAR IS COMMON TO PLEP AND DRAG. x 
HEHE HE HEH HEHE HEHEHE DE HE HE HE DE HE HE DEK HE DE HE HE HE HEHE DE HE HE HEH DE HE HE 2 DE DE HEE DE DE HE HE DE DE HE HE HE DE DE HE HE HE HE DE HE HE DE HE DE NE HE HE HEH EK 


5 
COMMON “7 SOLAR 7 F107 , FBAR , AKP , TVE 
INTEGER PSIZE, PSENS, PSYNS 
REAL MINSN 
x 
HE DE HEHE I HED HE HE DE HE DE HE HE DE HE DE HE 3 DE HE IE FEE IE HE DE DE HE HE DE HE IE DE HE DE DE HE DE FE HE DE DE DE HE HED DE HE DE HE HE DE HE DE 3E DE HE DE DE HEE HE EE HE HEH 
* INITIALIZING CONSTANTS. x 


HE HE HE HEH HEE EE EE HE EE EE HE EE EE EE EE EK ER HE EE HE EE EE EE REX 
DATA R2 / 541.15E-06 7 ,DEG /57.295779513D0 “ , ROE /6378.165D0/ 
See Ri Rng eT ae B eae nee Ga. ae 


* INPUTTING CONTROL PARAMETERS. 
1200500000002 2200000 2 0 II IID II III 


READ *, ITER1, ITER2, NOA, PSENS, IET, PSIZE 
HE HED HEHE HE IE HE 3E HE DE IE HEE I HE DEE HE 3 DE HE 3 DE HE HE 3 HE HE DE IE HE 3 DE HE 3 DE DEE DE IE HE DE DE IE DE DE DE IE HE DE DE IE DE DE HE DE DE DE HE DE DE HE HE DE DE HE HEE 


x INPUTTING SATELITE SPECIFIC PARAMETERS. x 
3 HD I I HEHE HEHEHE HE HEHE HEHEHE HE HEH BEI HE HE HE I HE DE DE 3 HE HE HEHE HEHE HE DE DE 3 3 3 3 DE 3 HE HE HEHE DE DE DE 3K 3 33K HE HE DE DE DE DE 3K DE 3 36 HE HE HE HE HE 


READ *, CD, A, D, SPI, W, WT 
x 
HEHEHE HEHE HEHEHE IE HE DEE HE HE HE DE HE DE DE 3 3 3K HE HE DE HE DE DE DE HE DE DE IEE DE HE DE DE HE HE HE HE HE HEE DE DE DE DE HE HE DE HE DE DE DE DE DE DE DE HE DEE HE HEE EEE 


X INPUTTING EARTH AND ATMOSPHERIC PARAMETERS. X 
HEHEHE HEHE K HEE HEHE HEHEHE HK HE HE HE HEHE HE HE HE DE DE DE HE HEHE HE DEK HE HE HE IE DE DE DE HE HE IE HEE DE DE HE HE DE HE DE DE DE IE HE HE HE DE DE DE DE DEE HE HEE 


READ * , F107 , FBAR , AKP , TVE 


PoP 2S PSS DHEL SS SES ELESES ESSE SE SESE SESE SES ES SSS ESE SE SESE SST STS SETS SSS FS SF 


- INITIALIZING BROUWER MEAN ORBITAL ELEMENTS FROM THE NAVSPASUR X 


x FIVE CARD FORMAT GRBITAL DATA SET. * 
1 HEHEHE HEH HEHE HE HEHE DE 5 DE HE DE DE SE FE DE HE FE DE DE SE 2 HEHE DEE 3 HE 3 3 HE 3 DE HE DE 3G HE DE HE DE IE DE DE DE HE DEK 3 DE BE DE DE DE DE HK DE DEK HEE X 
* 
READ 5 , UJD , ETU , HO , GO, ES, XI 
5 FORMAT € 78X,F14.8,5¢€1X,F8.49)) 
READ 10 , RND , ESD 
10 FORMAT (€ 20X, ake 7 LIX EL 2.) 
READ 15 , AO , ADOT 
15 FORMAT ( 78x, a 5,1X,E14.7) 


5 

HE DE HEHE HE HE HE DE IE HE DE HE HE HE HEE I HE DK 3 HE HE DE HED DE IE HE DE DE HE IK DEE 3 DE HE HE DE HE HE DE DE HE HE DE HE HE DE HE DE DE HE HE HEHE HE HE DE DE DE HE UE HEE HEE HE 
x INPUTTING THE REVOLUTION NUMBER AND CORRESPONDING CHANGE IN x 
*  SEMIMAJOR AXIS * 


HE HEHE HEHE HE HEHEHE DE DEE HE IE I SE HE HE SE HE HE IE HE HE EE EE SE DE DE DE EE DE DE DE SE DE DE DE DE DE DE SE DE EE DE DE EE EE HE EE KERR KKK KKK 


R 
TRC NOA .EQ@. 0) GO Gest 
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DO 20 K = 1, NOA 
READ *, IRVCX) , DACK) 
20 CONTINUE 


HEH HHH HEH HEH HEH HK HEH HE HEH HE HE HE HEHE HE HEH HEHE HE HK HEHE I HE HEH HH HEHE HK HE HEH HK HK HK HE HK HH HK HK KH HK 
* THIS NEXT SERIES OF CODE SETS THE INITIAL VALUES FOR THE SCANNING X 
% ARRAY MATRIX SO AS TO ALLOW THE RECOVERY OF THE ORIGINAL INPUT x 
a COCMei em Seno tT iVIiny MODE IS NOT SELECTED. X 
HH KH HE HEE HE EK HH HE HK HE HE HEH HE HE I HEHEHE HEHE HE HK HEHE HE HE HEHEHE HEH HK HE HEH HE EE HEH HK HH HK HK HK HHH HE KKH RK 


30 DO 90 ICN = 1,20 
SNSDELCICN) 1.0 
SNSMINCICN) 0.0 
SNSARCICN) = 0.0 

90 CONTINUE 


KKK RH KAR HK KK KH HK HH KKK HK HH KKH HHH HHH KH KKH HK HH HH KH HHH KKH KK HK HK KKH KK KKK KKK KKK 


* THE NEXT SERIES OF WRITE STATEMENTS OUTPUTS THE ORIGINAL INPUT * 
7 ~<) PARAMERERS. UNITS ARE SUPPLIED FOR CLARIFICATION. % 

HK HK HH HH HK HHH KKK AR KKH KKK HK KKK HH HH HH HK HH HH HHH HHH HK KK KH NHK AHH KK KKK RKKRKKR 

WRITEC6, 98) 

WRITEC6,99) 

WRITEC6,103) 

WRITEC6,99) ; 

WRITEC6,104) D, F107, FBAR, AKP, TVE 

WRITEC6,105) CD, A, SPI, W, WT 

WRITEC6,106) UJD, ETU, HO, GO, XI 

WRITEC6,107) RND, ES, ESD, AO, ADOT 

WRITEC6,108) ITER1, NOA 

IFCNOA .GT. 9) THEN 

WRITEC6,109) 


DO 27 NN = 1, NOA 
WRITEC6,110) IRVCNN), DACNN) 
2] CONTINUE 


END IF 
WRITEC6,99) 


HEH HE HE HE HEHE HEHE HE HK HEHE HE HEH I HE HE HE HE HE HE I HE I HE HE HE HE OE HE HE HE HE DE OE HE EE HE DE HE DE UE HE HE HE HEHE HE HEHE HEE EEK EK HK 
x THIS NEXT CODE SETS THE INTIAL INPUT PARAMETERS TO SOME RETAINER 
% VARIABLES THAT ALLOW RETURNING TO THE ORIGINAL VALUES IF THE BS 
* SENSITIVITY MODE IS SELECTED. THIS PREVENTS HAVING TO REREAD ALL * 
* THE INPUT PARAMETERS BACK IN FOR EACH RUN OF THE SENSITIVITY SCAN.X 
HEHE HH HEHE HEH HEHE HEHE HE HEHE HE HE HE HE HE I HEH HEHE HE HEH HEHE HE HE HE HE HE HEHE HE HE HE HE HEH HK HEH HE HEH HEH HH HK HHH HK 


V1=CD 
V2=A 
V3=D 
VG=SPI 
V5=W 
V6=WT 
V7=F197 
Vd=FBAR 
V9=AKP 
ViO=TVE 
V1Ii=ELJ9 
Vi2=cTU 
MIS =ne 
V14=50 
Vis =E5 
V16=XI 
V17=RUD 
V1I38=ESD 
V19=A0 
V20=ADOT 


HEH HEHE HEH HEHE HEE I HE HE DE HE HE HE HE BEE HE HE HE HE DE HE HE HE I HEH HE HE HE I HEH DE HE HE HEHE HE HEE HE HE HE HE HE HK HK KK 
* THIS IS THE ITERATIVE SCANNING ROUTINE. ARGUMENTS ARE INCREMENTED 
x FRACTIONALLY THROUGH THEIR SPECIFIED RANGE AS SPECIFIED IN THE x 
x Neue Trees onOUlINEMoNAL to <PoIZE LTTERATIONS OF PLEP>. 3 
Ae SENS 1 vee ODe £5 NOT SelECTED PSIZE=1. x 
HEH HE HEHEHE HH HEH HEH HEH HEHE HEHEHE HEH HH HE HHH HEH HE EH HE HEH HE HEH HE HE HE HEH HIE HE HE HEH HH HK HEH EE HEH HHA 
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TS aN -EQ. 0) THEN 


SIZE=0 
END IF 
DO 40 PSYNS = 0, PSIZE 
* 
PES SSESS ESSE SSS SSS SSS SESS ES ESS SES SESS STS SSeS eee es eee e sees ee eset ts tts ee: 
* THIS NEXT SET OF CODE IS ENVOKED WHEN PSENS IS EQUAL TO ONE. xX 


* IT SETS UP FOR AND CALLS THE SENSITIVITY ARGUMENT ARRAY SNALYS ak 
& SUBSEQUENT SCANNING. CPSENS = 1 SELECTS SENSITIVITY MODE 
59035333 ESO GEE I RI IEEE I DEBE 


TE *CPSENS ZEQ. 2 Saher 
TET=0 


x 

HEE HE HEE HEHE HEH HEHE I HE DE HE DE HE DE I HEE HE HEHE IE HE EH HE EH HE HEE HE HE HE HE HEH HE HE HE HE HE HE HK HE HE HEH HE HE HE HE HE HEH HHH KKK KK 

* THE SUBROUTINE SNALYS INTIALIZES THE SCANNED INPUT PARAMETERS TO X 

_% BE INCLUDED IN THE STUDY. I? ALSO PRENTS THewceio tty Ghahal ee cee 
HH HEHEHE HE HEHE I HE HEH HE DE HE HE HE DE I HE DE HE HEHE HE HE HE HK HE HE HE HE HEH HE HE HE HE HE HE HE HE HE HEE HE HEHE HE HE HE HE HE HEH HEHE HEH HE HE HK HH HH 


x 
TRCPSTNS CQ. OU mmimemn 
CALL SNALYSC PSIZE, SNSAR, SNSMIN, SNSDEL, SUMCK, DELSN, 
MINSN) 
ELSIE 
C 


D= Vi : 

HED HE HEHE HE HE HE DE DE HE HEE HE HE HE HED DE HE DE HE HE HE DE DE DE HE HE HEH DE HE HEHE HE HE HE HE HE HE DE HEE HK HE HE DE DE HE HE HE HE HE HE HEHE HE DE DE KK HEE HEN 
* RECOVERY OF THE ORIGINAL INPUT PARAMETERS AFTER THE FIRST RUN CF 
* A SENSITIVITY ANALYSIS. % 
HEHE HE HEHE DE HE HE HE DE HE HEHE DE HE HE DEDEDE DE HE HE HEHE HE HE DE HE HE HE HE HE HE DE HE DE HE HE HE HE DE DE HE HE DE DE HE HEE HE DE HE HE HE HE HE HE HE HEE EK IK HE EK KE 


A 


ZUM 
me] 
aa 


“Tl 

r—) 

Q 
FOu~s 4 uk et 
won 

< 

ae | 


FBA 


END IF 
END IF 
ASYNS=REALCPSYNS ) 
ASIZE=REALCPSIZE) 
IFCPSENS .&0. 1) THEN 
I> (SUMCK. EQ. 910) TBEt 
SCNFRC = ASYNSADELSN + MINSN 


SCNFRC = ASYNS/ASIZE 
ENO 
ENDS Fe 


HEH HE HEHE HEE DE HE HE DE HE HE DE HE HED HE HE HE HE HE HE HE DE HE HE I HEE HE HE HE HK HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HE DE HK HE HE HK HE HE HE HK HE HE EEK KE EK 
* THIS NEXT CODE CREATES THE CURRENT VALUE FOR THE SCANNED INPUTS. X 
x IF SENSITIVITY IS SELECTED THE PARAHETER WIUE SCAN -ROM SENG ae 
* WITH A GRANULARITY OF 1/PSIZE PER SCAN 

SS 3 EBD 3 ak NT 


CD SNSARC1) XCSNSMINC1)+ASYNS¥SNSDEL (1 


2 CDXCSNSAR(C1)-1.0 
A = PE Peete ih: 20 Sa ce SS 
= 4 


- 

Y-A ¥CSNSARC2)-1.0 
Y-D ¥CSNSARC3)-1.0 
J=SPEXCSHOAR CH) = 10 


SNSARC3) XCSNSMINCS)+ASYNS¥SNS DEL 


) 

) 
C3) 

SNSARCG) XCSNSMINCG)+ASYNS¥SNSDEL (4) 


) 
) 
) 
0) 
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Ae =m OAR CoD <C ONSHILNC 5) FAS YNSXSNSDELC5) )-W XC SNSARC5)-1.0) 

WT = SNSARC6)*CSNSMINC6) +ASYNSXSNSDEL (6) )-WTXCSNSAR(6)-1.0) 
FLO7= SNSARC7)*(SNSMINC7)+ASYNSXSNSDEL(7))-F107X(SNSAR(7)-1.0) 
FBAR= SNSARCE) XCSNSMINC38)+ASYNSXSNSDEL (8) )-FBARX¥(CSNSAR(8)~-1.0) 
AKP= SNSARC9) XCSNSMINC9) +ASYNSXSNSDELC9) )-AKPXCSNSAR(C9)-1.0) 
TVE=SNSARC10)¥CSNSMINC10)+ASYNSXSNSDEL C10) )-TVEXCSNSAR(10)-1.0) 
UJD=SNSARCI1J XC SNSMINCI11) +ASYNSXSNSDEL C11) )-UJDXCSNSARC11)-1.0) 
ETU=SNSARC 12) XCSNSMINC12) +ASYNSXSNSDEL (12) )-ETUXCSNSAR(12)-1.0) 
HO =SNSARC13) XOSHSMINC13) +ASYNSXSNSDELC13) )-HOXCSNSARC13)-1.0) 
GO =SNSARC19) XCSNSMINC14) +ASYNSXSNSDEL C14) )-GOXCSNSARC14)-1.0) 
ES =SNSARC15) ¥CSNSMINC15) +ASYNS¥SNSDEL(C15))-ES¥C(SNSARC15)-1.0) 
XI =SNSARC16)*CSNSMINC16 )+ASYNSXSNSDEL C16) )-XI¥CSNSARC16)-1.0) 
RND=SNSARC17) XCSNSMINC17 )+ASYNS¥SNS DEL C17) )-RNDX¥CSNSARC17)-1.0) 
ESD=SNSARC138) XCSNSMINC138)+ASYNS¥SNSDELC18)) -ESDXCSNSAR(18)-1.0) 
AQ =SNSARC19) XCSNSMINC19)+ASYNS¥SNSDELC19) )-AOXCSNSARC19)-1.0) 


ADOT=SNSAR( 20) ¥(SNSMIN(20)+ASYNS¥SNSDEL (20) ) -ADOTXCSNSAR(20)-1.0) 


ES SS TT eee Te ee eee ee ee ee SS STS SSS SSE SSS SSSSSSSS SSS SESS Sess sesssesssessesss 
* DEL IS THE CONSTANT MULTIPLIER FOR THE INTEGRAL IN EQUATION 335 OF x 
* NSWC TR 83-243 * 
TESS SSESSSSSCSSSSSESSSS SS SS SSS SS SSS SSS SSS SS seSsessssessesesessessesssss esses! 


DEL = CDXAX3.14159D0/SPI . 
x 
PES ETELSESLS SESE SESS SS SSS ESE SS SS SESS SSS SES SSS SS SESS EES SESE SES ESS SES ESS S: 
<p oN ovr ANALYSIS MODE ITS NOT SELECTED, THEN THE SYSTEM 
* PRINTS STANDARD MISSION LIFE TIME TABULAR DATA. x 
HEH HEHE HHH HE HEH HHH HEH HEH HE HE HH HEHE HE HEHE HE HE HK HEHE HK HEH HE HE HE HE HE HE HK HEH HE HE HE HE HE HE HE HE HE HE HE HH HK HH HK KH HK HK 


x 
Deer SENS EG.) 0) THEN 
BCO 


F =CD XA 
WRITE (6,98) 
WRITE (6,99) 
WRITE (6,100) 
WRITE €6,99) 
WRITE (6,101) BCOF, SPI, W 
WRITE (6,102) F107, FBAR, TVE 
ENDOLE 
HH HH EH HEH HEH HH HEH HHH HEHEHE HE HEH HE HHH HE HEH HEH HH HHH HEH HE HEH HEHE HE HEHE HEH HE HEH HEHE HEH HH HK HEH HK HEH HK HK 
* CONVERTING DEGREES TO RADIANS x 
HH HH HEH HEH HHH HEH HEH HHH HEHE HHH HH HEH HE HEH HEHE HE HE HE HEHE HE HEHEHE HK HEH HEH HE HE HEH DHE HE HE HE HE HH HE HEHE HEH HEH HEHE 
x 
XI = XI 7 DEG 
GO = GQ ” DEG 
HO = HO / DEG 


FLO = ETU / DEG 
* 


PE SSS SSS eS SOPeSreSeESesSEPSSSSSESESSS SSS SESS SESE SSE SSE ESS S SS SESS SESE SES SS SS 
elatomniex) CODE SEQUENCE CONDITIONS TRE ELEMENTS AS PER NSWC TR 82-3387 
* PAGE 7. IT CONVERTS THE KAULA SEMIMAJOR AXIS AO IN EARTH RADII 10 x 
* BROUWER MEAN SEMIMAJOR AXIS VIA EQNS 2.1 AND 2.2. THE SEMIMAJOR AXIS 
* DECAY RATE ADOT IS CONVERTED IN A LIKE MANNER VIA 2.5, 2.4%, AND 2.5 X 
HEHE HE HEE HE HE HH HEH HEH HEHE HEH HEH HEH HEHEHE NE HK HL HEHE HE HE HEHE HE HEN HE HK HE HEHEHE HEHE HEHE EH HE HEHE HH HE HEHE EE HE HK HK HK 
x 


FN =1.27¢1.- ES ® ES )¥¥1.5D0 
os =esrto( 2x) 2 
TB = TA * TA 
fe = le = €S.4IB) 7 2. 
TAD = (€€ 3.84R2) 7 €2.*A0¥A0D) )X¥TEXFN 
AK = AQ 
AQ = CA0XCC1. + 2.*TAD) 7% C1. - TAD) )¥*¥0.6666666/7D0)*ROE 
XDOT = TAD¥(3.X*ESDXCES/(€C1.-ESXES)) - 2.*CADOTZAK) ) 
ADT = CADOTZAK) XA0 + 2.¥AKXROEXXDOT¥C CL. + 2.*TAD)X® 
al CCl. - TAD) ¥® 5))¥*¥(€-0.3333333333D0) 
ESS SES SSS SSS SES See ee Te STS SSeS ee ee STS SESS SST eee eee sss esses ssessesssseseess 
* TIME IS NOW CONVERTED FROM SECONDS TO DAYS % 
PT SSS SSS SS ST ESS SSeS SSS TS TTS eee SS SS SSSI SSS SSS SSeS SSS SSS SSS SS SS SSS SSS eS ee: 
ADT ADT 7 86400.D0 


ESD ESD 7“ 86400.D0 


137 


RND = RND¥2.D0%3.14159D0%C1.7086400.D0)**2) 
eas 


HEHE HE HE HE HEHE HE HEHE I HE HE HE HK HE HE DE DE DE DE DE HE DE DE DE HE HE HE HE HE HE HEHE HK I HE HEH HK HH HEHE HK HH HHH HHH HH HEH HE HEHE HK EE HK HK 
* THE SUBROUTINE KOZAI IS THE OPERATING WORKHORSE FOR THE PROPELLENT x 


* LONGEVITY COMPUTATIONS. IT ALSO DECREMENTS PROPELLENT AS IT IS USED 
HEHE HE HEH HK HE HHH HEHE HEHE HEHEHE HEHEHE HEHE HE HEHE HEH HE HE HE HE HEHE I HE HEH HE HE HE HEHE HEHE HEHE IEE HE HE HE HE HE HEH UE EE ESE HK 


x 
CALL KOZAICADT,AO,ES,XI,FLO,G0,HO,UJD,W,DEL,ITERL,ITER2, 
SCNFRC, D) 
% 
HE DE HE HE DE DE DEE DE HE HE HEHE I HE HE HE HEHEHE HEHE DE HE HEH HEHEHE HE HE HEH HE HE HE HE HE HE HE I HE HEH HH HK HH HK HE HEH HE HEHE EEK EK HE HK 
* FORMAT STATEMENTS USED IN THE PLEP INITIALIZATION PROGRAM. % 


HHH HK AHH HK HH HK HK HE EEE EK KE HHH EK EK HEHE KEK KK HK KEE EK HEE HK EK EEE EEE EE EK 
X 


98 FORMAT ( 1H1 ) 
99 FORMAT(//5X, '!------------------------------------- ----------- Mee 

100 FORMAT(//15X, "PROPELLANT LIFE ESTIMATOR',//) 

101 FORMAT(//5X, 'CDA=',E12.5,' SQ KM'/5X, "SPECIFIC IMPULSE=', 
XF7.2,' SEC'/5X,'INITIAL WEIGHT OF PROPELLANT=', 
XF7.2,' KG!) 

102 FORMAT(//5X,'F107=', F8.2/5X, 'FBAR=',F8.2/5X, 
X ' TIME OF VERNAL EQUINOX=',E22.14) 

103 FORMATC15X, ‘INITIAL INPUT CONDITIONS') 

104 FORMAT(///1X, "ATMOSPHERIC INPUT PARAMETERS'//5X,'D (RELATIVE ATMOS 
XPHERIC ROTATION) = ',F10.46/5X,'F107 (DECIMETRIC SOLAR FLUX) = ', 
XE7 .2/5X,'FBAR (90 DAY AVG. F107) = ',F7.2/5X,'AKP (GEOMAGNETIC IND 
XEX) = ',F6.3/5X,"'TVE (TIME OF VERNAL EQUINOX PASSAGE) = ',F16.9, 
X' MODIFIED JULIAN DAYS'//7) 

105 FORMATC1X, "SATELLITE PHYSICAL PARAMETERS'//5X,'CD (DRAG COEFFICIEN 
XT) = ',F5.2," UNITLESS'/5X,'A CCROSSECTIONAL AREA) = ',F11.9,' SQU 
XARE KILOMETERS'/5X,'SPI (MOTOR SPECIFIC IMPULSE) = ',F8.2,' SECOND 
XS'75X,'W CINITIAL FUEL MASS) 2 ',F9.2,' KILOGRAMS'/5X,'WT CINITIAL 
X SATELLITE MASS) = ',F11.2,' KILOGRAMS'///) 

106 FORMATC 1X,"NAVSPASUR ORBITAL ELEMENTS'//5X,'UJD CEPOCH) = ',F14.8 
X,' MEAN JULIAN DAYS'/5X,'ETU (MEAN ANOMALY) = ',F8.4,' DEGREES'/5X 
X,'HO (MEAN RIGHT ASCENSION OF THE ASCENDING NODE) = ',F8.4,' DEGRE 
XES'/5X,'GO (MEAN ARGUMENT OF PERIGEE) = ',F8.4,' DEGREES'/5X, 

X'XI (MEAN INCLINATION) = ',F8.4,' DEGREES') 

107 FORMATC5X,'RND CRATE OF CHANGE OF MEAN MOTION) = ',F11.9,!' REVOLUT 
XTIONS PER DAY SQUARED'/ 7 
X5X,'ES (MEAN ECCENTRICITY) = ',F8.4,' UNITLESS'/5X,'ESD (ECCENTRIC 
XITY DECAY RATE) = ',F16.9," 1/SECONDS'/5X,'AOQ CKAULA SEMI-MAJOR AX 
XIS) = ',F11.5,' MEAN EARTH RADII'/5X,'ADOT (SEMI-MAJOR AXIS DECAY 
XRATE) = ',F16.9,' '///) 

108 FORMATC1X,'CONTROL INPUT PARAMETERS'//5X, "MAXIMUM REVOLUTIONS PER 
XITERATION *,17,/5X, "NUMBER OF SCHEDULED ORBIT ADJUSTS TO SEMI-MAJO 
XR AXIS) Se 

109 FORMAT(//5X,'SEMI-MAJOR AXIS INCREASE SCHEDULE'//5X,'REV ADJUST’, 
X5X, "CHANGE IN KILOMETERS'Z) 

110 FORMAT(5X,17,7X,F6.2) 


40 CONTINUE 
= OP 
END 


2 22 2222 22.222. 2,22, 2,2,2,2,2, 2,2, ¢, 2,2, 5,. 5,9, 9) 
HHH HEHE EE EE EH EK EE ERE KKK KKK KKK 


SUBROUTINE SHALYS %X 


SOO IOI OUI ODIO 
e222. 2.2. 2.2.2.0, 2.2.2.2, 0, 2,0, 5,2, 0,2, 2, 0, 2,0, 5,0, 5,5,1 


mK OK OK OOK OK OK OK OOK OK OK KOK OK 
*K 
mm 


Serre age ae SNALYS CPSIZE,SNSAR, SNSMIN, SNSDEL, SUMCK, DELSN, 
MINSN) 
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KK OK OK OK OOK OK OK OK 


DIMENSION SNSAR(20), SNSMINC20), SNSMAX(20), SNSDEL(20) 
CHARACTER¥G SNSCAR 

LiVEGER MIE PSIZE 

REAL SNSAR, SNSMIN, SNSMAX, SNSDEL, SUMCK, DELSN, MINSN 
WRITE (6,98) 

WRITEC6,120) 

WRITEC6,109) 

WRITEC6,120) 

WRITEG6G, 200) PSIZE 

ASIZE=REALCPSIZE)D 

READC 5,102) 

Suc. = 0.70 

DO 20 N = 1, 20 


READ 101, SNSCAR, SNSARCN), SNSMINON), SNSMAXCN) 
SNSDELON) = CSNSMAXCN)-SNSMINON)IZASIZE 
IF (SNSARCN) .EQ. 1.0) THEN 
SUMCK SUMCK + 1.0 
DELSN SNSDEL CN) 
MINSN = SNSMINCN) 
WRITEC6,300) SNSCAR, SNSMINCN), SNSMAXCN), 
l SNSDEL CN) 
END IF 
20 CONTINUE i 
Te CSUueK Gh. .1.0) THEN 
DELSN = 
MINSN 
ENDe LE 
WRITEC6,120) 
WRITEC6,400) 


98 FORMAT ( LHL ) 
120 FORMAT(//5X, '------------------------------------------------ 127) 
100 FORMAT(//17X, "SENSITIVITY ANALYSIS FOR',/16X,"LOW EARTH ORBIT ', 

X "SATELLITES',/19X, "PROPELLENT LONGEVITY',//) 


200 FORMAT (5X, * TABLE OF INPUT PARAMETERS SCANNED FROM MINIMUM TO’, 
X 75X,' MAXIMUM WITH THE INDICATED INCREMENT PER RUN. ',75X, 
X * THE GRANULARITY IS 17',14,'.'/75X, * ALL OTHER INPUT PARAMETERS 
X ARE RETURNED TO THEIR',/75X,* INITIAL VALUES AT EACH INCREMENT. '! 
X “4//*INPUT PARAMETER!',8X,*MINIMUM',13X, "MAXIMUM',13X, ‘INCREMENT 'Z) 


102 FORMAT (1X) 

DOLeGORMA © ASameLX, FS.15 0616.9, .9X, F16,9) 

300 FORMAT (5X, AG, 6X, 303%, F16.9)) 

G00 FORMAT(//7, 10X, ' SENSITIVITY SCAN RESULTS FOR ABOVE SPECIFICS'// 
15X,*REV NUMBER',8X,*TIMECDAYS)',15X, "FRACTION OF SCAN‘) 


RETURN 
END 


oor 


0. 
0. 


PES SEES ELSES SESS ESSE SS ES SES ESS | 
HEH EH HE HEHE HHH HE HEH HEHEHE HEH HE HEHE HE EEK 
xX «x 
XX SUBROUTINE KOZAT BS 
xx x 
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a 


KOK OOK OK OK OK OK AR OK KOK OOK OK OK OOK OK OOK OOK OOK OOK OK OOK OK 


HK HHH KKK KKK KK KKK KK KKK KK KKK KKK KN 
PES SSS SESS SESE SESE SSS SESS SESS Si 


PS SSeS SSE SESE SES ESSEC ESSE SESE SC SESE SESS SESS SSS SESS SS SSS SS SESE S SSS SSE S SSS 3 
*® 
SUBROUTINE KOZAI IS THE SYSTEM WORKHOUSE. IT CALLS BRAUER TO CON-x 
VERT THt BROUWER MEAN ORBITAL ELEMENTS OF THE NAVSPASUR INPUT TO 
BROUWER OSCULATING ELEMENTS. THESE OSCULATING ELEMENTS ARE CON- 
VERTED TO KOZAI-LIKE MEAN ELEMENTS IN THE SUBROUTINE MEAN. THE 
SUBROUTINE PERIOD COMPUTES THE INITIAL ANOMALISTIC PERIOD. THE 
SUBROUTINE GEOP COMPUTES THE EFFECTS OF THE OBMATE Geol) Jie ine 
ORBITAL ELEMENTS CTHROUGH THE J4 TERM). SUBROUTINE THRST COMPUTES 
THE EFFECTS ON THE ELEMENTS DUE TO PERIGEE BURNS TO CHANGE 
THE SEMIMAJOR AXIS DURING IRVCK) BY THE AMOUNT DACK). FINALLY, 
SUBROUTINE DRAG COMPUTES THE THE AMOUNT OF FUEL NEEDED TO OVERCOME 
THE DRAG EFFECTS BY CONTINUOUS INTRACK. MICROTHRUSTING. 
KOZAI KEEPS TRACK OF THE FUEL USED AND ORBITAL ELEMENT CHANGES. 
WHEN THE FUEL IS EXHAUSTED KOZAI RETURNS SYSTEM CONTROL TO PLEP. 


CHE HEH HK RK HH HHH HH HHH HHA HHA HANK HK KK HK HK KKK KKK KKK KKK K KK KKK KKK KK KKK KKK 


KOK OOK OOK OOK OOK OOK OOK OK OOK OK OOK OK OK 


SUBROUTINE KOZAICADT,A0,ES,XI,FLO,G0,HO,UJD,W,DEL,ITER1,ITER2, 
a SCNF 


RC, D) 
DIMENSION XOSCC6) , XMC6), DXDTC6) , XDDC6), IRVC1O), DAC10) 
DIMENSION DXiMC6) 
COMMON 7 OADJ/Z NOA, IRV, DA, SPI, WT, PSENS 
DATA BO,B2,B3,B4,B5,B6 7 .398603254E+06,-.1755528999E+11, 
iL -26386647738E+12, .10630735996E+16, .805605022E+18, 
1 -(/2921158 5665047 
DATA DEG/57. 295779513D0/,CN2/0.0/7,A1,E1,RN170.0,0.0,0. aie 07 
DATA PI 7 3.14159D0 / 
DATA NPLOT 7 0h 


IRC = 0 
PCOUNT=0 
KK = 0 


Pile =.2.00x% ia 
CALL BRAUER (B0,B2,B3,B4,B5,DT,A0,ES,XI,FLO,G0,HO,CN2,A,CE,CI,CL, 
1 G,H,A1,E1,RN1 ) 


IF (PSENS.EQ. 0) THEN 
WRITEC6,100) ITER] 


100 FORMAT(//5X,'MAXIMUM NUMBER OF ORBITAL REVOLUTIONS=',I110) | 
WRITEC6,101) ITER? 
Tei FORMAT(//5X,'DATA PRINTED EVERY ',I3," ORBITAL REVOLUTIONS') 
SCNFRC = 0.0 
END IF 
CALL PERIOD € AO, ES, XI, PA ) 
XOSCC1l) = A 
XOSCC2) = CE 
XOSCC5) = Cl 
XOSCC4) = G 
AGSCCS)) =H 
XOSCC6) = CL 
CALL MEAN CXOSC,XM) 
ILINE = O 


CIW = XMC3)*DEG 
GN = XMC4)%DEG 
HW = XMC5) DEG 
CLW = XMC6)XDEG 
IF CPSENS SEG. 0) SIREN 
WRITEC6,102) 
ro eon CXMCK),K=1,2),CIW,GW,HW,CLW 
END 
103 FORMATC/5X,'SEMIMAJOR AXIS=',E22.14,' KM'Z 
X 5X, 'ECCENTRICI 1) =Se2cn1 47 
X 5X, 'INCLINATICN=!,E22.14,' DEG*'Z 
X 5X,*ARGUMENT OF PERIGEE=',E22.14,' DEG'/ 
X 5X,'ASCENDING NODE=',E22.14,' DEG'/ 
X 5X,*MEAN ANOMALY=',E22.14,' DEG') 
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x 


1o2 FORMAT(//5X, "KOZAI MEAN EPEMENT SET") 
CALL THRST ¢ oe DXM, IRC, KK, PA, DWOA) 
DO 10 i= OER 
\ s2@ 1 <= Y ) % DELT 7 86400. 
TT = UJD + T 
he CALL GEQP (XM, DXDT) 
16 DO 20 JJ = 1,6 
XMCJJ) = XMCJJ) + DXDTCJJIXDELT + DXMC JJ) 
Geers . Li. 4 ) GC TO 20 
AMC JJ) = AMOD € XMCJJ), PIe ) 
20 CONTINUE 


DO 30 II = 1, NOA 
Hee EQ eR YCET)) IRC = I 
Peer Re .cQs 1) KK = If 
30 CONTINUE 
CALL THRSTCAM, DXM, IRC, KK, PA, DWOQA)D 


ott 2.2.2.2 2.0 2.0.2, 0,2,2,5,2,9,5,2,0,2,5,0,%,0, 2,9 %,0, 0, 0,0, 0,5,5,5,5,5, 9,5, 0,0, 9,5, 4,0,0, 4,990, 9)5,5,0. 55,8 tetetitetatatetetatates 


* CALCULATING THE DEL1 IN NSWC TR 83-243 EQN 11. % 
HEHE HN HEH HH HE TENE HEHE HE HE HE IE HEH HE HEHE HEH HHH HE HEH HEHEHE HH HEH HE HH HEH HEH HK HK HK KH HHH HK HK K 


¥% 


KOK MK OK MK OK 


RP=XMC1)*¥C1.0-XMC2)) 
VP=SQRTCCBO/XMC1))*CC1 .0+XMC2) 9701 .0-XMC2)))) 
DEL1=(1.0-CRP/VP)*DXB6XCOSCKXMC 3) ) ) X*2 
DDELT=DELXDEL1 


CALL DRAG ( XM, DDELT, TT, DELMG ) 
W = W + DELMG - DNOA 
WT = WT + DELMG - DWOA 
ERCP OENS) 68.50) THEN 
IF (PCOUNT.NE.CITER2-1)) pach 
PCOUNT=PCOUNT#1 


: GO TO 10 
PCOUNT=0 
ENDL 


Deere eine GE S500) ifLiINE = 0 
Tete = ett NeG +] 
TE € ILINESSEQR IO OSWRITE € 6,104) 
Peo [lPNEPeeGy 2) ARi ee. 6,105) 
WRITE € 6, 106) I, TT, T, W 

a, LG a 

ME Ce e000) THEN 
DEO CNRE ORS. 69. LIEN 

PRINT ®, T, SCNFRC 


WRITEC6,107) I, T, SCNFRC 
END IF 
SeUURIN 
END IF 
10 CONTINUE 
FE SCHPLOT: 69. 1) °1HEN 
PRINT *, T, SCNFRC 


E 
WRITEC6,107) I,T,SCNFRC 

END IF 
104 FORMATC1H1) 
105 FORMAT(/5X,'REV NUMBER',3X,'TIMECMJD)',8X,'TIMECDAYS)', 8X, 

X ‘PROPELLANT REMAININGCKG)',77) 

DUG RORHIATCSA, [10,3X%,F9.2,46xX,F10.2,15X,F8.2) 
107 FORMATC5X,110, 8X,F10.2,15X,F16.9) 

RETURN 

END 


pete eett reese trsitetshetttests. 
KHKHKK KKK NAKA KK HK KR RKKRKKKKK KER 


xX %X 
*X SUBROUTINE BRAUER XX 
4 XX 
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KK KK 


HK HE HHH HHH HH HH HK HK HK HK HE KE EH 


SUBROUTINE BRAUERCBO,B2,B3,34,B35,DT,A2P,E2P,CI2P,CLO2P,G02P,HO02P, 
1CN2,A,CE,CI,CL,G,H,ADT,RND, ESD) 


1004 A2P2=A2Pxx2 


AZ2PG=A2P2X*2 

CNO= SQRTCBO/CA2P2¥AZP DD 

C2ZP2Z=EZPH¥Z 

ETA=SQRTC1.-E2P2) 

SINEI=SINCCI2P) 

THETA=COSCCI2P) 

THETAZ=THETAX¥2 

THETAG=THETAZ#X2 

THETA6=THETAGXTHETA2 

CJ2=-B2/02. ¥BOXA2P2) 

ETAZ=ETAX¥2 

ETA3=ETA2ZX%ETA 

ETAG=ETAZ%¥2 

CJ21P=CJ2/ ETAG 

CJUSIP=83/ CBOXA2P2ZXACPXETAGXETAZ2) 
CJGIP=(3.¥5B4)7(08.XBOXAZPGXETAGXETAG) 

CUSIP=B5/ CBOXA2ZPGXA2PKETAG*X2*XETA2Z) 

FUN1L=3.*THETA2-1. 

FUN2=1.-5.x*THETA2 

SINEI2=SINEI¥*¥2 

AL=A2P¥CJ2*FUNI 

AQ=-A1L/ETA3 

A2=3.*A2PXCIJZ¥SINEI2 

FUN5=1.-11.*THETA2-(€G40.*XTHETAS)/FUN2 

FUN6= -~FUNI~-C8 .XTHETASIJ/FUN2? 

FUNG=THETA2/SINEI2 

FUN22=FUNZX*2 

CJZIP2=CJZ21P**2 

EQIP=(C CE2PXETA2) *¢3. ¥CJ21P2*FUNS- TO. 
ICJGIP¥FUN6 37029. *CJ21P) 

EZle-=—-2 .xEUPP 
TNT S banana cies naienaaeee) Tie gO EES SNE I) 

J2iP) 

El1LP=-.75D0*XESIP+¢CC. 25D0XETAZCESINET ) XC GUSI1P+ Si Zag0xCJ51iP 
1*(€G4.+3.x*E2P2)x% 
2(¢1.-9 .*THETA2-( 24. *THETASG)/FUN2Z) D)/CJZIP 
CIO=-CE2PXIHETA)/CETAZ +o NEL)D 

CIZ2=CJZIPXTHETAXSINEI*¥1.5D0 

CII=E2PXCI2%.66666667D0 

FUN7=C-.5DOXETA3*¥CJ21P)/E2P 
CLZIP=CETA3/CJ2Z1P)%C . 25D0XCJZIPZ¥FUNS— .833353533535D0XCJGIPXFUNG ) 
CLI2P=CNOXC1.+1.S5DOXCJZIPXETAXFUNI+.09375D0XCJZIP2ZXETAXC~-15.4+16. 
L¥ETA+25. XETA2+030 .-96. XETA-90. XETAZ) XTHETAZ2+°(0105.4+1449¢. *XETA+25. 
2¥XETAZ) XTHETAG) + .9375D0XCJGIPXETAXE2P2*(3.-30.*THETA2+35 . ¥THETAG) ) 
CL22P=.5D0XCNOXCN2 
G21P=C€1./€294.4CJ2Z1P) 2¥C-3.*CJZIP24C2. tE2P2-11.*02.4+3.xXE2P2)¥*T KHETAZ 
1-G0.%€2.4+5.*E2P2) XTHETAG/FUNZ-G0O0. XECPOXTHETA6/FUN2Z2)4+10.*CJG1IPX 
2Gee 
34+£2P2-3.¥C€ 2.45. E2P2) *THETA2-8.*(2.4+5.XE2P2) XTHETAG/FUN2-80. XE2P2*% 
GTHETA6/FUN22)) 

Ay Py ena aiaieieaeitbe Maria 2 SD EN AES 

xET 2 
20690.-192.¥*ETA-126. XETAZ) XTHETA2+0385.+360,. *XETAt+G5. XETA2Z) XTHETAG) + 
3.3125D0XCJGIPXC21 .-9 . XETAZ+0 270.4126. XETAZ) XTHETAZ4+06385.-189. 
GRETAZ) XTHETASG)) 

H2=1.5DO0*CJZIPXTHETA 

H3=-2. HZ 

H1=.66666667DOXE2PXH2 
ASIP=(€C35.*CJSIPXE2P2XE2PXTHETAI/C1494.*CJ21P)I¥¢€ .SDO/SINEITXCFUN?-C 
ar Rie aie oe ras oN 0A ISON 
HLIP=-.25% 
] H31P+¢€.25D0XE2ZPXTHETAI/CCIJZIPXSINEL) )XCCJ31P+.3125D0*XCJ5SI1PX 
PEW reseed tiara ae ea ISN 2 LIN 207) ote 

XSINEI 2x 
46(04.+3.¥*E2P2)*(3.4+€16.*THETAZI/FUN2+(40. XTHETAGI/FUN2Z2) ) 
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0510 


Bau 
530 


100 


200 


210 
220 


8018 


8022 
8023 


ne lP=CECRZXIHETA) A Glia XC J ZIP) ¥(-3.xCJZ1P2*C11.+030. XTHETAZ)/FUN2+ 
SET Mee er nae C16 AT MEDACIZPUNZEEGO. RTHETAS D7 
eS ideeg rae ee” URS TOD PERC 5 DETREI ATS -RETACY 

eee 

236. XETA-5. XETA2) XTHETA2) +1. 25D0XCJG1PXC5.-3.XETA2Z) *(3.-7. XTHETAZ2) ) 
AID=CJ51P/CJ21P 

AID2=FUN2-(16 .XTHETAG)/FUNZ2 

C1=35.7384. XAIDXETASXE2PRSINEITXAID2 

AIDS=THETA2/SiNEI 

AIDG=THETA2XSINEI 

Ege seecrexeer 

C2 selec. xe LDC C-EZPRSINEIXCS +2. XE2P2)+EZP3XAIDS)XAIDST 
12.¥*E2P3¥AIDGXC5 +032. XTHETAZ)/FUNZ+(080 . XTHETAG)/FUN2Z2)) 
C35=1.7-9.*THETA2Z-€29.%THETAG)/FUNZ 

AID5=CJ31P/CJ21P 

C4$=. 25D0*AID5X%(-E2PXAID3)+5.764. XAIDKC-E2PXAID3%(49.+3.XE2P2)+ 
LE2PXSINEI ¥(26 .+9. XE2P2) )¥C3-15./32. XAIDKXEZPRAIDGRC4 43. XE2P2)% 
203.4016. THETAZI/FUN2Z+090. XTHETAG )/FUNZ2) 

Co-beh7 cl telAs) xCs =e cP 2x03 .-E2P2)) 

C6 -Ceer ease +61 XGe2eexE2PZ 07 CCG Foe eecP ZTE ARCS. +9 XE2P2) >} 
C7=.25DO0XAIDSXSINETXC5+5 464. *C3XAIDXETACHSINEI C6 
eee re eatin UNMIS RS TNE EACH TD REE CIRCE 

=D 

CL2P=CLI2PXDT+CL22PXDTXX2+CLOZP + RNDXDTXDT 
CL2P=AMODC(CL2P,6.2831853071796D0) 

IF(CL2P)520,530,530 

CE2P=CUZr+6.283518550/1796D0 

G2P=G12PXDT+G02P 

H2P=H12PXDT+HO2P 

SINEG=SINCGZP ) 

COSING=COS(G2P ) 

DIUE=STHEGACSINEGRCESTPXOTNEGSEZIEOSELIPIFEDIP 
HIP=CCH3IPXSINEG+H2Z1P)XSINEG+H11P) XCOSING+H2P 

GEBRoecr seber yo UUACCE CIP +G2lP)XSINC 2. xGZ2P)+(C1+C2)*COS(3.% G2P) 
1+(C4+C7 )¥COSING 

CEIP=CL2P 

U=CL2P 

DELTAU=CU-E2PXSINCU)-CL2P)/ C1. -E2PXCOSCU) ) 

U=U-DELTAU 

IFCABS(DELTAU)-1.€-10)200,100,100 

ae ni CUE2PRSINCU)-CL2P)/C1.-E2PCOSCU) ) 

SLHETR=SINCE) 

COSE1LP=COSCE) 

GIR -G2P 

ADIVR=1.7(01.-E2P*COSEIP) 

SINFIP=ADIVRXETAXSINEIP 

COSFIP=ADIVRXCCOSELP-E2P) 

FIP=ATAN2ZCSINFIP,COSFIP) 
DECABSCFIP-CL2P)-3.1415926555e9ab0)220,210,210 

STOP 

FUN3=C1.+CN2*T)*®.6666666/7D0 

COUSEG=COS(2. XC(GIP+FIP?) 

Slwmo-ornee.<(GIPtFI1P)) 

ADIVRZ=ADIVR*XZ 

ADIVR3=A4DIVRKX3S 
Gl-Cl2RrOro-UVETCIIXSINFIPXSINEGS( 2. xXClIxCOSFIPtCI2)xCOSrFG 
FUNSSFIP-CLIP+E2PXSINFIP 

H=H1P+€2.8H1LXCOSFIP+H2) XSINFG-HIXSINFIPXCOSFG+HSXFUNS 

KFUN=H/6 .28513853071796D0 

FUNO=KFUN 

H=H-FUN9X6 .2831853071796D0 

IFCH)8022, 8025, 3023 

H=H+6 .2831853071796D0 

A=A2P/FUN3S+A0+CAL+AZXCOSFG)XADIVR + ADTXDT 
AID6=ADIVR2XETAZt+ADIVR 

AIDT=SINCZ. XGZPTELIE) 

ATDS-SINCZ.XG2ZP+3. xF1P) 

D1=.2500XCJ2Z1PRC6 .XC5 . XTHETAY—1. -)XFUN8+(3. -5.*THETAZ)*C3. XSINFG+ 
13.*E2PXAID7 +E2PXAID3 )) 
D2=.25D0XCJZ1P*C2. 4035. XTHETA2Z—1 . DXCAID6+1. )XSINFIP+3. XC1.-THETAZ)# 
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KKK KOK KOK OK OK MK OK 


KK KK KK OK OK OK OK OK 


1 CCmWALIDGO41 . DRAILD/ FCAITUG+ . 33533555500) vs 
AIDS=COSCZXGZP TRE. 

AIDIO=SCOS CZ axGZPr a XEITPS 

D35=-ETAZ®. 5DCXCJZIPXC1 .-THETAZ) *C3.xXAIDI+AIDI0O?D 
ETA61I=1./CETA3SXETAS) 

DS4=ETAGIXCCS® +COSTIPXCS . tECEXCOSR IE XCS tect -COSR TE. 
DE=ETAGIXC c2P+COSFIPR(3 .+E2PXCOSFIPXG3 3 ecee ens ope 

D6= ETAZXCJ2%.5D0%(C35. XT HETAZ=1 .) XD4435 x Cl eA ADS COSEGeens 
GAL=GPLP+D1+( E2PXETAZ)7 CL eA 

CE=CE2ZP—1. XC 1. +CN2ZXDT D¥X .66666666666667 DUT] ORE UG tee ona 
EDL=.5D0XE2ZPXCLZIPXSIN €2.%G2Z?)+C3ExXCOS INGTEZCPACI Mes Se Gze | 
lL ETAS*D2 

AIDIS=SINCCL2P) 

AID15=COSCCL 2P } 

ESL=CE *AID14+EDLXAIDI5 

ECL=CE XAIDIS5-EDLXAIDIS 

CEP Ont ECl XECU+ ESEXESL D 

CL=ATAN2CESL,ECL) 

G=GAL-CL 

G=AHODCG,6.2831853071796D0) 

IF(G)3024,3025,8025 

$0294 G = G + 6.28351855077 I60e 
8025 Suen 
ND 


PESESELESSSES TS SSS SSS SSS SSS FFF 
PSSST SESE SSS SES ESS SESE SESE SS SF 
XX HX 
HX SUBROUTINE POSVEL %X 

XX 
HEHE HE HE HEHE HE HEHE HEHE HE HE HE HIE IE IE HE HE HE HEHEHE HH HEHEHE 
HEHE HE HE HE HE HE HE HE HE HE HE HE HE HE HE HEH HEH HE HE HHH HH HEHE 


SUBROUTINE POSVEL CA,CE,CI,CL,G,H,BO0,R1,R2,R5,V1,V2,V5,R, VEL) 
CALL NWTRPH (CL,CE,E,SINEP,COSEP) 
HSIN=SINCH) 

HCOS=COSCH) 

GSIN=SINCG) 

GCOS=COS(G) 

CISIN=SINCCI) 

CICUS-CUsKer 
ALL=HCOSXGCOS-HSINXCICOSXGSIN 
AlZ2=-HCOSXGSIN-HSINXCICOS*®GCOS 
AZL=HSINXGCOS+HCOSXCICOSXGSIN 
A22=HCOSXCICOS¥GCOS-HSINXGSIN 
ASL=CISINXGSIN 

AS2=CISINXGCOS 

FUN=SQRTC1.-CEXCE) 
RL=AXCALIXCCOSEP-CEJ+ALZ*CFUNXSINEP DD 
RZ2=AXCAZ1XCCOSEP-CE)+A22*XC FUNXSINEP DD 
R3=AXCAZ1LXCCOSEP-CEJ+A32XC FUNXSINEP DD 
PaAx( (CE SCocen. 

FUN LSCSQRTCBROXA ISR 
VISFUNIXCALILXC-SINEP)+A12XFUN*COSEP 3D 
V2=FUNLACAZLT XC -SINEP D+A22XFUNXCOSCP ) 
VS=FUNLECASL RC -STHEP JtA32XFUNXCOSEP )D 
VELSFUNLTESQRTC1. -CCEXCOSEP ) X¥2) 


RETURN 

END 
PES SESS SESESTESE SESS SS STE SESS SS: 
HE HE HEHEHE HE HEH HEH HHH HH HEHEHE HEH HE IE IE HE HE HC HE HE HE 
XX XX 
KX SUBROUTINE NWTRPH XK - 
%X *X 


PES Sesser ses essere tse sersss ese: 
HHH HHH HH KH EHH HHH HK AHHH HH KKK KKK 
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OK OK OK KOK OOK OK OK OK OK OOK OK OK OK 


KK KOK OK OK OK OK OOK OK OK 


SUBROUTINE NWTRPH (€QL,B,E2,SE,CE) 

NEWTON RAPHSON SOLUTION TO KEFLERS EQUATION SUBROUTINE 
EL=QL+C(CB¥SINC QL) )701.-Bx¥COS(QL)) 

SE=SINCEl) 

CE=COSCE1) 

E2=cl+CQL+BXSE-E1)/(C1.-B¥CE) 

IFCABSCE2- ElJ= 1.E-8 )4,4,3 


SE-StN< EZ) 
CE=COSCE2) 
RETURN 

END 


HK HH HH HK HHH RH HHH HHH HHH HH MH KKK 
HH HK HHH HH HK HH HH HH HH HK HK HHH HK NRHN 


XX XX 
HX FUNCTION XSPA HX 
xX KX 


HK HHH HHH HK HHH HK HK HK KK KK RK MK KR 
HHH HK KK HH HK HK KKK HK KK KK KKK KKK KAR 


THE FOLLOWING FUNCTIONS DETERMINE THE FIRST ORDER SHORT 
PERIODIC VARIATIONS AS A FUNCTION OF THE MEAN ELEMENTS. 


FUNCTION XSPACXM) 

DIMENSION XM(6) 

REAL JC4) 

DATA J,R71.,0.00108276D0,-0.00000255D0,-0.00000156D0,6378 .165D0/ 

F = XMC6)+02.*XMC2)-XMC 2) X¥3..7G4+5./96 .XXMC2) XS. DHSINCXMC6 DD) 

A ete eae cm ace na cco ee th! 7292. AXMCZI AKG. )XSINCZ. X 

XM 
Meno + (135.712. xXMC ZIRK2. 45.7192. UC? Faas) JESINGS. XXMC6)) 

P = XMCI)XC1 . -XMC 2) XZ.) 

Ramee P7CL xNC2)xXCOSCF)) 

XSPA= JO2)*RER“XMCLIXCOCXMCLIZ RR HHS. 4001.73.72. XSINCXMCO3) SIN 
LOXMC3)) 43.72. XSINCXMC3) )XSINCXMC3) COS C2. XCXMOG)+4+F) ))- 
201.-3.72.*¥SINCOXMC 3) DRXSINCXMC3)))*C1. -XMO2) X¥2.) RKC -1 .5D0)) 

af 


HH HHH HH HHH HK HHH HHH HK RH HK KKK KK KKK 
HK KH KH KKH HK HK HH KKH KKK HK HK HK KKK 


% x XX 
¥X FUNCTION XSPE HX 
XX ¥%X 


HH KH KKK KKK HHH HARK KH KH KA KH KKK KKK 
HK HK HHH HK HH HK KH HK RK HK HK KKK KK 


FUNCTION XSPECXI4) 

DIMENSION X?tHC6) 

REAL JC4) 

DATA J,R71.,0.00108276D0,-0.00000255D0,-0.00000156D0,6378.165D0/ 

eS XMCSIFC2]. XXNC2) -XNC2) HHS. 744+5./96 . XXNC 2) KS. DHSINCXMCG) ) 
Perc ODOR NIG exe 2, = 110724. XXMC ZENG. +17 7192. *XMC2)¥RG6  JRSING 2. X 
2xM¢ 
36) )+C13.712. XXMC 2) XX2.-G3.7192. XXMC2)*X5. XSING 3. XXMC6)) 

P = XMC1) ¥C1.=XAMC 2) 8X2. } 

XSPE = JO2)/2. CRSP )XR2. KC]. - 1. SDOXSINCXMC3) DESINCXMC3) )) 
Bie ee aang ences SAIC 2) XRT SDOI/AMCZI4S. XC 1. + 
2xMC2) 
SRXINC2974.)*COSCFI+1 . SDOXXMC 2) XCOSC 2. XF+XMC2)*¥XMC2)79.*COS(3.¥F)) 
443.785 . XJ C2) €CR/P ) RZ. KSINCXMC 3) DESINCXMC 3) *CC1 411.79. XXM02) X¥2. ) 
5 ¥COSC2. ¥XMCGI+FX4XMC2) XXMC 20749. XCOS C2. XXMCG)-F +5. *€XMC 2) ® 
GeO Sean 4 TEXZ C7 ty oa RANK 2) EXMC2))73.*COSC 2 . ¥XMC9)4+3.XF) 
7+1.5D0*XMNC2) *COSC 2. XXANCG) 44. HF +XMC 2) XXMC 2074. XCOS C2. ¥XMC4)4+5.*F) 
8+1.5D0%Xi4C2) XCOSCS. XXMCG))) 
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KK KK MK MK OK MK OK OK OK 


KK KK MK OK OK OK OK OK OK 


WK OK OK OK OK OK OK OK OK 


RETURN 
END 


HHH HH HH HH HHH HHH HE HEHEHE HH HE EE HEE HHH 
Pe Se SSS eSESSSLSESE SE SESE SES SE SE SF | 


xX xx 
HX FUNCTION XSPI HH 
XX XxX 


PP Pe Set SPS SSSE SESE SSE SES SESS FF F 
PSS SSSSESESSESSSE SSE SESS SS SSE SSS SF 3 


FUNCTION XSPICXM) 
DIMENSION XMC6) 
REAL J(4) 
DATA J,R41.,0.00108276D0,-0.00000255D0,-0.00000156D0,6378.165D0/ 
F = XMC6)40€2.%*XMC2)-XMC2) 3.7445 .796 . *XMC 2) *5. DKHSINCXMCG )) 
LECT 25 DORAN C 2c ci i 29 a gee 
36))+€13.712. ¥XMC2) ¥¥2.-43.7192 .¥XMC2)**5. SINGS. XMC6)) 
P = XMC1)¥*C1.-XMC2) 2.) 
XSPI = 3.78.*502) *CR/P) 2. HSINCZ.XXMC3) I XCXMC2)¥COS( 2. ¥XM(4)+F) 
L +COSC2.%CXMCG)+F) )4XMC2973.*COSC2.*XM(4)+3.%F)) 
ges 
ND 


HHH HHH HHH HH HH HH HH HHH HK KK KK KK KK 
KKH KKK KKK HK RK KRHA RK KKK KAKAKRAKKKAKK 


xX xx 
xX FUNCTION XSPW xX 
XX XX 


HH HH HH HH KK HK KH RK KKK KK KKK KR KKKKK 
HH HH HH HHA HH KH KH RHR KKK KRAKRKRKKKKKR 


FUNCTION XSPWCXM) 

DIMENSION XM(6) 

REAL J(G) | 

DATA J,R/41.,0.00108276D0,-0.00000255D0,-0.00000156D0,6378 .165D0/ 
F = XMC6)402.*XMO2)-XMC 2) 43.4445 .796 . RXMC2) HHS. DRASINCXMCE DD 
ed Ae a tae cea Cilia gig i SANS 
2XM 
269)40135.712. XXMC2) RK2.-G3.4192. XXMC 2) HHS. DRHSINGS. XXMC6)) 

P = XMC1)¥C1.-XMC2) ¥¥2.) 

XSPN = 3.74. 4502) €CR/P) RR2 HOG. -5. ASINCXMC3) DXSINCXMO 3) )) 
LCE ACG FAME ORS EN SP Rigiger OSAP a amet hoa 
RAL 

i Soa aaa aaa a a Sea NS 
4/12 

5.)-1. 5D0¥IC2) CRSP) HR2 RC CSINCXMC 3) )XSINCXMO3) 7G. XMC 2) RRZ 2. RO 


61.-15.78 .XSINCXMC3) ) 2.) /XMC2) RSINC 2. RXMCGIFFIFXMO2)716 SIN 
7 CXMCS) 442. RSINCZ. XXMOG) —F I+. SBOXC1. 2. SDOeSTHCKIICS ID 442 XS INC 


8 CXMCG)+F)2=-€7.712. XS INCXMC 3) 0 ¥¥2. -XMC 2) X10 2976. XO 197s. XSIN 
9CKMC3) 2 %¥2. VIZXMO2]IRSING 2]. XXMCG)43.4F)—-3.78. XSINCXMC3) R22. XSIN 
1102. €XiAC4 +4. KF MXC 20716. XSINCXMCS) ) 442. XSINC2]. RXMCG)45 . KF) D 
12-9 .716 . XJ CZ) CRSP RR]. RSINCXMC 3) ) RR2. RHSINCZ. RXMCG) D 

Saas 

ND 


HK KK KKK KKK KKK RK KK KKK KK KKKKKKKN 
HE HH HH HEH HHH HEH HEH HH HHH RH KKH HH R KK 


%X x 
xx FUNCTION XSPO XX 
%X %X 


HEHE HHH HHH HK HHH HE HE HEHE HE HE HE HH RR 
HH HHH KH TEKH KH HK HRA KKK KH KKK KKK KKK 
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2K OOK OOK OOK OOK OOK OOK OOK OK OOK OK 


2K OK OK KOK TK OOK OOK OOK OK OK 


FUNCTION XSPOCXM) 

DIMENSION XM(6) 

REAL JC49) 

DATA J,R41.,0.00108276D0,-0.00000255D0,-0.00000156D0,6378 .165D0/ 
F = XtMC6)+02.4XMC2-XMC2) ¥¥3..7945 796 . XXMC 2) X¥5. DXSINCXMCG )) 
By aoe emma eo 82 F727 192 AKMCZ)EXG . )¥SINCZ .X 
36))40135.712.XXMC2) X¥2.-G43.7192 . XXMO2) XH5. DXSINCGS . XXMCE DD 

Parte 1 XC 1. RKC ZAR?) 

XSPO=-1. 5D0¥JC2) XCR/P)*R2. XCOSCXMC 3) XC F-XMC6 I +XMCZIXSINCE) -XMC2) 
1/2. ¥SINC2].XXMCGI+FI-SINC2. XOXMOGI4F2 72. -XM02)76 .XSINC2. XXMOG4)43.% 


HHH RHR HH HHH HHH AHH HK RKRKAKXRAKKRN 
HHH HH HHH HARRAH KK KH HR KKKAKARKKKRKK 


xX xx 
+s FUNCTION XSPM HX 
xX PS 3 


HRA HARA RAK AR RARAARAKARAARAKARAARRARRR 
RHKRRRRAKARRKK KR RKRARARRARRRARRKARARR 


FUNCTION XSPMCXM) 
DIMENSION XM(C&) 
REAL J(4) 
DATA J,R41.,0.00108276D0,-0.00000255D0,-0.00000156D0,6373 .165D0/ 
F = XMC6)402.*XMC2) -XM0C 2) %*3 2.7945 .796 . XXMC2I)*R5. XSINCXMCG) 2D 
1+€1.25D0*XMC2)**2.-11 4264. XXM02) HKG. 417 07192. XXMC 2) RG. IRSINCZ.R 
2XMC 
36 33*C13.712.XXMC2) X¥2.7-93.7192 .XXMC2) XK¥5. DASINCS . HXMC6)) 
P = XMC1)4C1. -XMC2Z)24%2.) 
aay = -1.5D0XJC2)*CR/P) 42. *SQRT CL. -XMC2) 4¥2. 7 XMO2)4CC1 = 1. SDOX 
SINC 
Se KGimmae ys tn soe I CORBSENC EI EANC CI 72 ROTM CAPO 
XMC2 
ee a eS GS) ASTNCXMC SDORCC~  ODGIKC1 «4 
1.2500 3 
6%XMC2Z)*XMC2Z) DXSINC2. XXMCG4I4F 2 -XMC2) XMC2)78 .XSINC2.XXMOG)—-FD47 276. 
nme Sa AE) TODORAMC 2 ESENCE AAMCS IES. 
8x 
9I+XMC2) 42.73 .KSINC2Z. XXMOG4)45.%F)I)+9.716. XS CZIRCR/P DV RR2.X 
ISQRT C1. -Xi402) ¥¥2. DRSINCXMC3S) DRSINCXMCS) DXSING2. XXMCG) 
es ) 
ND 


HAH RKA RA RRARKARRAARRARARARARARAARR 
HH MRK HHH HHH MH RH HMA ARKARRARKRARKARR 


XX XX 
xX FUNCTION XSPR xX 
4X x 


MRR MAH KK RAK AKARKRK RR AK KK KARA RRARKAKR 


FUNCTION XSPRCXM) 

DIMENSION XM(6) 

REAL JC4) 

DATA J,R41.,0.00108276D0,-0.00000255D0,-0.00000156D0,6578.165D0/ 
F = XMC6)402.*¥XMCZ2)-XMC2) XS .7945 796. XXMC2) HRS. DASINCXMCE DD 
1+€1.25D0*XMC2)*¥2.-11 2.724. XXMC2) KRG. 417 7192. XXMC2) HKG. DASINGZ.X% 
2XMCG6)34013.712. XKMO2) ¥2.-G435.7192 . XXMC2) RK5. DASINGS.XXMCE) 

P = XMC1)401.-XMC2)*X2.) 

RRe= P7CL.tAhC2)xCOSCE) 

oe - .SDOXJC2Z) XC R¥¥2.7PI*XC1.-1.5D0XSINCXMC 3) XSINCXMC3)) 0401. + 
27€1.+SQRTC1.-XMC2)%82. 2 XCOSCF 42.7 CSQRTC1. ~XMC2) 42.2 DXRRZXMCL 
3) I+. 25K%IS C2) XRHX2./PRSINCXMC3) D2. XCOSC 2. XXMOG4D 42. FD 
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mK KK MK OK OK OK OK OM OK 


KK OK MK OK 


KR KK KK KOK OK OK MK 


10 


“20 


RETURN 
END 


HK KH HK HHH HK HK HHH HK HH KKK KKK KKK KKK 
HK HH RH KAR H HK HHH HHH NHK KK KHK KKK KKK 


%X XX 
€X SUBROUTINE MEAN x 
XX xX 


KK HHH KH HHH HK HHH KH HK KKK KKK RK KKK 
HK HHH HM HHH HHH HHH HHH KK KKH HK HK MK 


SUBROUTINE MEANCXOSC, XM) 


COMPUTES THE MEAN ORBITAL ELEMENTS FROM THE OSCILATORY ELEMENTS 
AND THE FIRST ORDER PERIODIC VARIATIONS BY THE RELATIONS, 
XM = XOSC -XSPCXM). 


DIMENSION XOSCC6),XMC6),XMNC6),XSP(C6),TOLC6), DELXMC6) 

POUT es amine cic erica aaecaal ee Un 
KOUNT = A 

INITIAL GUESS FOR MEAN ELEMENTA. 

DO 10 1=1,6 


XMCI)D = XOSCCTI) 

CALCULATE THE SHORT PERIODIC VARIATIONS. 
XSPC1) = XSPACXM)- 

XSPC2) = XSPECXM) 

XSPC3) = XSPICXM) 

XSPC4) = XSPWCXM) 

XSPC5)3 = XSPOCXM) 

XSPC6) = XSPMCXM) 

KOUNT = KOUNT+1 

DO 50 J=1,6 


XMNCJ) = XOSCCJI-XSPCJI) 


DELXMCJ) = ABSCXMNCJI-XMCJ)) 

IFCDELXMCJ) .GT.TOLCJI2 XMCJ) = XMNCJ) 

CONTINUE 

DO 60 K=1,6 

IFCDELXMCK) .GT.TOLCK) .AND-KOUNT.LT.30) GO TO 20 

CONTINUE 

IFCKOUNT.GT.30) WRITEC6,100) DELXM 

CONTINUE 

FORMATC/T5, 'CONVERGENCE OF ONE OR MORE MEAN ELEMENTS WAS NOT MET.' 


I7'THE ABSOLUTE DIFFERENCE IN THE ELEMENTS ARE SHOWN. '//7E12.5) 


RETURN 
END 


HH KK HRHKAK KAKA KKK KKK KK KK KAKRKKRKKKR 
HK HHH RHR KK KKK KKK KKK KKK KKKKKKN 


xX XX 
RX SUBROUTINE GEOP %X 
xX xX 


2E 5 HEE KE HE EK EK HE HEHE HK EN EK HK HK AAR RKKHEX 
HK MH HK HHH RK KKH HHH AK AKHK KAR KKAARK 


SUBROUTINE GEOP CXM, DXDT) 

DIMENSION XfMiC6) , DXDT($6) 

REAL N,NCONS, NCONST, NCON 

REAL JC4) 

DATA J,R/41.,0.00108276D0,-0.00000255D0,-0.00000156D0,6378.165D0/ 
DATA U/398600./ 

N = SQRTCU/CXMC 1) XX3.9) 

P = XMC1)¥C1. -XMC2) 882.) 

DXDTC1) = 0. 

S3 = SINCXMC3) DXSINCXMC3)) 

l SAIC 2) Xxc- 

COSCXMC3)) 

DED = C-3.9732.4NKJC2) S02) CRSP) RRG RS SRC14. 7-15. S35) XXMC 23% 


>< 
“ 
a 80 


IXS¥SINC2. XXMC4) ) 
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KOK OK OOK OOK OOK OK OK OK AK OK 


DEDO = 3.78. ¥*NKJC 3) XCR/PI ARS. XSINCXMC3) ¥0G.-5.S3)*XS¥COSCXMC4) D 

DEDOT = 15.732. ¥*NXJCG) RCR/P) XRG. XS3RCG6 M7. XS3)KCXMO2)¥XS) 
LXSINC2].¥XMCG)) 

DXDTC2) = DED-DEDO-DEDOT 

DID = 3.464.¥NXJC2)RSC 2) CRSP) ERG. RSINCXMC 3) 2.) *014.-15.%S3)% 

UX 2) 42. 4SINC2]. ¥Xi104)) 
Ao > 3.78. RNKICSIRCR/P) RRS. ACOSCXMO3))*04.-5.%S3) XXMC2)*COSCXMCS 

DIDOT = 15.764. ¥NXISCG)RCR/P) ERG. RSINC 2]. XXMC3))*C6.-7.¥S3) XXMC2) 
IXMC2]) *SINC2Z2.XXMCG)) 

DXDTC3) = DID+DIDO+DIDOT 

DND = .75D0XNXJC2)*CR/P)X¥2.404.-5.%S3) 

Oa 6 OG = 1, 4535) 455-028 . = 158 45341355. 4S 342. IAXMC 2) XZ. 

DNDO = 3.716 .¥NXJSC2)¥42. CRSP) RG. H0G8-103.*S3 + 215.74. 8S 3KK2 2+ 
Be ete DU Rosato Ocoee eA 2) eee ro. XC]. —1, SD0xSS) 404. - 5.53) 
2 SQRT CXS) -.25D9*XDW¥COS (C2. ¥XM(C4))) 

PRINT %, XMC2), XMC3), XMCG),R,P,N,J03),583 

DADGT = 3.73. 4NAICSIECR7PIAASEC(S .-5 . XS3 IRC SS-XMOC 2) XX 2ZHCOS (XMS) 

A Rn bi ame an amma ery °F 2b. FS SIEAMC2) DAS ING AM 

Du = SSAC. =e =O neo, KOO tol. SRSSRXC. )RXMO 2) RXZ. 

POO) l= 52 ie Ct OR AKG KOLO. Oe. XS 5T99 .XS35RK2Z.+.75D0%024.- 
Tote So STO Ss. tOoe ke . VEX 2IARS. tDWNXCOSCZ. XXMC4))) 

DXDTC4) = DWD+DNDO+DWDCT-DWDCTT 

DROS =C=1- SD0UANAIC 2) XCR/P RAZ ORCS 

DODO = 1.5DOXNXJC2)¥%2.*CR/P) ERG. RCSKHC2. 250041. 5¥SQRTCXS)-S3 
DACA UU Hee SDUXSQRIC XS) EXMC 2) RAZ .79.8%61.41.25D0X53)-XNC2) 2X2. 
27 AX 7. -15. 4535) *COS C2. XXMC4) )) | 
Hates). 5.7 Sm ( SIACRAP XA SR AC 15. XS5-4 RXMCZ)¥CS/SINCXMC 3) DESING 

XMC4) 

DODOTT=15./16 . XNXJ CG) RCR/SP DERG. RCSHC C4. 7. HS3IKC1L. $1. SDORXMC 2) 8X2. 
Cato. at. Xo5 ) RXMGZ) REZ. HEOSCKNC41¥2.)) 

DXDTC5) = DOD-DODO-DODOT+DODOTT 

NCON = 3.78. 4NXIC2) 8X2. 

NCONS = NXJ0G)¥CR/P DERG. 

NEGNST “=535.738 .ANXIJCSIJACR/P) RS. 

DM = N¥C1.41.5D0*JC2) ¥CR/P)€¥2. 401-1. SDOXSS)RSQRTCXS)) 

Ne mee ee ere 0 as 31 DBRS 5 1A 5AR Ret 1 22R0xC 1. 
2.5D0% 

Teele Oko noo) eon Ai ey ome SX 5) XXMCZ)¥*24535/16.%C14.- 

Z Do- ea 504 Cl 2 oUUAMIIC 2) X42) XCOS C2. XXMC4)) DXSQRT CXS) D 

Da (See U Aon 4) oe KOOKS e Ol ODOTD .XSSt1L1L? ./6 . RaS¥S35IRXMCZ)RRZ 
eat ee aoa OL. 7S. S558 55 175, XXMC2)ERG. ) 

Wee = AtlGe anc .76.X95%(70.-125,455tC 996 .-06 .¥S53IRXMC2) RZ KCOS C2. 
LXi1049) 427.7128 . RXI40 2) HKG. RSSRSSRCOS CS. KXMCG) ) 

DMDO = NCONKACR/SP)¥RG.SCSQRT CXS) ¥C3. X¥DM1+DM2) 

PMDGT = HOONSTASENCXMNCS) 0404. -5.XS3)%C1.-9.. XXMC 2) REZ. D/AMC2)® 
ISQRT CXS) ® 
IS INCXM094))-465.7128 . RNCONS*C3 .-G0 .¥S34+35.XS3RS3)¥XXMC2) R¥2. XSQRT CXS) 


DMDOTT=15.764%.NCONSKS5¥(6 .-7. XS3)C2.-5. XXMC 2) 4X2. XSQRT CXS) COSC 


12. *XMC4)) 

DXD+(6) = DM+DMD+DMDC-DMDOT+DMDOTT 
Sas 

Esl 


KAKKKKKK KKK KKK KKAKKKKKKKKKRKKERR 
RKRKKRKKK AK KRRAK RK RA RKRKRKKRKXKER 


*X xx 
RX SUBROUTINE SOLOC XX 
KX xX 


KKK KK RK RAK KKK ARKKKARKKKKKRKRARKR 
KRRRKRKK KKK KKK KKK RK KKK KKKKKKKKR 


SUBROUTINE SOLOC ( TVE 


fees, OLA > DS » AS » AB ) 
Peas crei 25.500 7 , PI'7 


3.14159D07, TOL’ 0.0001D0 / 


Pil =s0— SDOXPI 
Be, = 2. 8 PI 
pas = 15004 PI 
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AR OK OK OOK OR OK OK OK OK OK OK mK OOK OK OK OK OOK KOK OK OK 


mK OK 


RPD 
DEL 
ERS 
CLO 
TE 

IFC 
a5 
DS 

BLA 


DO 10 | ee 
J=z=tTe-1 
ARG = 0.5D0 * J xX PI 
TST = ARG = CLON 
Totl. =<ABSC 15a) 
DESC TSil. « LE 4) Iie Col Geco 
10 CONTINUE 
ETA1L = TANCDS) 7 TANCEPS) 
ETA = ASINCETAL) 
ETA = ABSCETAD 
Te €C CLON-. LT . PIL) « ANDSex eLUlia wc 
IFCCCLON.LT.PI).AND.CCLON.GT.PI1)) AS = PI - ETA 
IF (CCLON.LT.PI3).AND.CCLON.GT.PI)) AS = PI + ETA 
IF CCCLON.LT.PI2Z) ANDECCLON SGI Fis eAs eee ieee 
GO TO 21 
20 AS = CLON 
21 A5 = AS + BLAR 
RETURN 
END 
PSS SSSS Ese ses sess sesesssesesess: 
pees eee ese sess esses sesse sss rs sss: 
%X *xX 
¥X SUBROUTINE DATPRP % x 
%X XX 
SSS SSS SSeS sssessessssesesesesst se: 
PSST STS SSE Se Sess sssessesescsessesss: 
‘SUBROUTINE DATPRP (€ X,Y,XX,YY ) 
COMMON 7“ PREP 7 ASUN, DSUN, RAS, DECS 
X = ASUN 
Y = DSUN 
XX = RAS 
Woe=nOGECS 
RETURN 
END 
pee TESTS SSS SESE Ss Sess esseseseess 
eee Se Sess ses ese sess sess esses se: 
XxX *X 
4X SUBROUTINE PERIOD XX 
%X XX 
Pe SSS TSS TESS SESS SS OS SHES Ss Sees: 
HK HHH HEH HH HK HE K KKK EK KAN 
SUBROUTINE PERIOD CoAG SES, xl 
REAL wc 
COMMON 4“ XXXX 7 IET 
DATA J2 “ 0.00108276D0/, RZ 6378.165D0/ , U4398600.7,PI73.14159D0/7 
PA = 2.¥PIXCSQRTCAOQXA0DKA0/U) )¥XC1.203. XJ 24RER¥C3.KXSINCXIDVESINGXI 
] ~2.0)/04. XAODXA0*(C1.-ES*®ES)*¥*1.5D0)))) 
IF (€ IET . EQ. 0 DJ RETUR 
WRITE (6,100) PA 
100 FORMATC/10X,"INITIAL ANOMALISTIC PERIOD= 
BeTURS 
ND 


= PI / 180.0D0 
T=T - IVE 
= EPS1 * RPD 


N = 0.985647D0 * RPD * DELT 

CCLON.LT. 0.0 ) CLON = Pie + Chan 

CLON.GT.PI2) CLON = AMODCCLON,PI2) 
= SINCEPS) * SIN (€ CLON ) 

= ASINeGesDSs > 

R = BLA ® RPD 


HK HK KH HHH KKH HHH KK HK HK KK RA KRKKRE 
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if | 


020)) AS a= ETA 


SEC", 77) 


KK KOK OK OK OK OK 


KK KK OK OK OK OK OK OK OK 


HE HK HHH HEHE HE HE HEHE HK HEE HE HE HE DE EE HEE EEE KEK 


xx %*X 
*% SUBROUTINE SALT XX 
xX xx 


HEHE HE HH HE FE HEX HEHE HE HEHE HE EE HEE EEE EEK HK EH 
HE HK HE HE HEH HEHE HEHEHE EE HEHE KE KK HEH HEHE K RK HK 


SUBROUTINE SALT (D,XM,RM,Z,RE ) 

DIMENSION XMC6) 

DATA ROE/$373.165D0/, EQE2/6 .69342E-03/ 

MemomK Oe sohncl. +CEQE2/( 1. -E0EZ)JxSINCD)XSINCD)) 
Ceeenhitet Aor RAM) — RE 

RETURN 

END 


HH HEH HK HH HK HK HHH KH HK KKH KKK KKK KKKKR 
PES eS SSeS TESTES SESE SL ESS SST ESS FS 


XX xX 
£% SUBROUTINE SATLOC xX 
XX xX 


HK KH HK HK HHH HK HHH HK HK KKK KK KKK KKK 
HE HE HE DE HK HEE IE EE EE DEE EE EE EE EE KEKE MME 


SUBROUTINE SATLOC (MTYPE, XM, DELT,DSUN,ASUN, DECS,RAS,GMGLT,R) 
DIMENSION XMC6) 

DATA BO/ 398605. 254D0/,W.7292115E-04/, PI/3.14159D0/ 

PIO2 = 0.5D0xPI 


Rie = 2.*P1 

P2 = PI2 

GeOo = o7eb 72. 

AG = AMOD € WXDELT , Pl2 ) 

A = XM(1) 

CE = XM(2) 

CI = XMC3) 

G = AMOD € XMC4) , PI2 ) 

H = AMOD € XMC5) , PI2 ) 
CALL POSVEL(CA,CE,CI,AM,G,H,B0,R1,R2,R3,V1,V2,V35,R,V) 
DECS = ATAN.C Bos an + R2*R2 )) 


RAS = ATAN2 € R2, RI 
feet MirPemeea 00) GOLTOC1O 


AEP = TANCOSUN)D * COSCCI) 7 SINCCI) 
ALP = ASINCALP) 
ALP = ABSCALP) 
DECS = DSUN 
ee CoM PE EQ. ye DeECS = = DSUN 
Ul = ASIN SINCDSUN) 7 SINCXMC3))) 
U2 = ASIN € SINCDECS)/SINCXMC3))) 
Tl = ASUN + PIQ2 
T2 = ASUN - PIOZ2 
Wil = G 
OMA = H 
Per (t2 .LT. G20*) GO 10.20 
=<.Tl .6T.. PiceeGOmiOes 0 
Te COOMA LE. Tig eAnoe ew OMA «SE. 1T2))"GO TO 40 
GO TQ 50 
20 elo = PIlZ + i2 
IF (CCOMA.LEVTID CAND. COMA.GE.0.0)).0R.(C(OMA.LE.PZ).AND.COMA.GE.TS 
l 2d) GO 
GO TO 50 


SOUS oe a, - Pi2 
DECC COMMAS Getic) - al COMA, PE, Pil2))-OR.(COMA.GE.0.).AND.COMA.LE.TS 
l 2.) GO TO 40 
GO TF sq@eCla GE ._PI0Z ) GOTO 41 
TEC) DECSaet |) . SCC. ey RAS = OMA - ALP 


41 RAS = OMA - ALP 


Loa 


MK WK OK ON OK OOK OK OOK OK OK OK 


KK OK OK mK OK OK 


* OK 


IF € DECS . Lie 0. Oar RAS =s0M4 + one 

IF € MIYPE . E@ . 1 ). RAS = JRAGe + rE 

RAS = AMOD ¢€ RAS , PI2 ) 

GO TO 10 

IF € CI . GE . PIO2 3 GO TO 5l 

RAS = OMA + PI - ALP 

IF ¢DECS . LT . 0.0) RAS = OMA + PI + ALP 
GO TO 42 

RAS = OMA + PI + ALP 

IF ¢€ DECS . LT . 0.0 ) RAS = OMA + PI - ALP 
GO TO 42 

ELNG = RAS - AG 

IF CELNG .LT. 0.0 ) ELNG = PI2 + ELNG 

ELNG = AMODCELNG,PI2) 

GMGLT = ASIN € 0.9792D0*¥SINC DECS)+0.2028D0*COS( DECS)*COSCELNG - 
] (291.. XPI71302300 

RETURN 

END 


HH HRA KKK KKK K KKK RR KKRKAKK KKRKKRKAKKK 
ttt. 2,2.¢,. 2,0, 0,5,9,4,0,5,5,5,5,0,0,0,5,5,5,5,5,5,5, 5,53 


xX xX 
XX SUBROUTINE DRAG %X 
xX %X 


HR RAR KH RK RK RK RRAR RARKRKARKRARKRKKR 
Lt +t. 2,5,0,4,0,,5,0,4,5,,5,0,5,5,5,5,5,5,5,0)5,5,5,5,5, 2 


SUBROUTINE DRAGC XM, DEL, TT, DELMG 3) 
REAL N,1I0,11,12,13,19,15,16,17 
COMMON 7 SOLAR 7 F107,FBAR,AKP,TVE 
COMMON 7 PREP 7 ASUN,DSUN,RAS,DECS 
DIMENSION XMC6), SUNCG) , SATC2) , GEOC3) , BFC8) 
DATA BO0/398603.254D0/,BLA/30./,PI/3.14159D0/, EOE’ .003355D0/ 
DATA IYR41981D0/7,TVD/1.7,THR/17 .7,TM414.7,T7S/56.24¢3D0/ 
I 


Pig e526. Aer 

SUNCS). = Za. Pata ) 
SUNCG4) = 1. 

GEO(1) = Fi07 

GEOQC2) = FBAR 

Ged¢3) = AKP 

N = SQRTCBO/CXMC1) ¥3)) 
.TCOR = AMOD € XMC6) , PI2 3} 
TCOR = TCOR / N 

Wy TT - € TCOR 7 86400. ) 


“TTR = 1721044. + 367XIYR - (7%IYR“4) + TVD + (CC THR¥3600. + TMX60. 
1 + TS ) 7 86400. ) - 2400001.0D0 

DELT = Til ser 

stot bell tan Iaith aS NOG): 

MTYPE = 0 

CALL SATLOCC(MTYPE, XM, DELT, DSUN,ASUN, DECS,RAS, GMGLT,RM) 
STR = XMC6) 


X'4€6) = 0.00 
CALL SALTC DECS,XM,RM,HP,RE ) 
XMC6) = STR 


CALL DATPRPC SUNC1),SUNC2),SATC1),SATC2) 2) 

CALL ATMDENCTTT, SUN, SAT, GMGLT, HP, GEO, RH0, DRDH, SH) 

RHOP = RHO * € 10. ¥X¥9. ) 

SoP = SH * 0.001D0 

SOypeoeny Me » DSUN,RAS,ASUN, BLA,HP,F107,RHOP,SHP) 
MTYPE = 

CALL SATLOCCMTYPE, Xi, DELT, DSUN,ASUN, DECS, RAS, GMGLT,RM) 
CALL DATPRP © SUNC1),SUNC2),SATCO1),SATC2) ) 

CALL ATMDENCTTT,SUN,SAT,GMGLT,HP,GEO, RHO, DRDH, SH) 
RHOMIN = RHO * € 10.%%9 

SHMAX = SH * 0.001D0 

OUBELa oS” aan » DSUN,RAS,ASUN, BLA, HP, F107, RHOMIN, SHMAX) 
M = 2 

CALL SATLOCCMTYPE, XM, DELT, DSUN,ASUN, DECS, RAS, GMGLT,RM) 
CALL DATPRP € SUNC1), SUNC2), SATC1), SATC2) ) 

CALL ATMDENCTTT,SUN,SAT,GMGLT,HP,GE0,RHO,DRDH,SH ) 
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RHOMAX=RHO * €10.%*9 ) 
SHMIN = SH*¥0.001D0 
CALL JAC60(DECS, DSUN,RAS,ASUN, BLA, HP, F107,RHOMAX, SHMIN) 
A = SINCDSUN)XSINCXMC3) )XSINCXMCG) )+COSCDECS )*¥CCOSCXMC5)-ABLG)X 
1 COSCXMCG) 3 -COSCXM(3)) XSINCXMC5)-ABLG) XSINCXM(G) )) 
B = SINCDSUN) XSINCXMC3) )¥COSCXMC4) )-COS(DECS) ¥( COSC XMC5) -ABLG) 
1 SINCXMCG4))+COSCXNC3) )XSINCXMCS )-ABLG) XCOSCXMC4G) )) 
F = CRHOMAX - RHOMIN) 7 CRHOMAX + RHOMIN ) 
RHOO = RHOP - 0.5D0XCRHOMAX + RHOMIN) * F XA 
= XMC1) % XMC2) 7 SHP 
emo) = € FRA 7 C1, + F X A D) 
Ceram fF xa) —~ 2% € 1. = 0.5D0KZ) = 0.500 
Wie SZ Area 2X 2KZKZZ 
CTHAVG = 0.0 
IF € TM1 .GT. 0.00 ) CTHAVG=(Z1+SQRTCTM1))7(€2%Z) 
SHNM = ABSC SHMAX - SHMIN ) 


N 
_ 
bb us 


HAVG = 0.5D0XCSHMAX+SHMIN)+ € SHNM 7 CSHMAX+SHMIN)) XAXCTHAVG 

BETA = 1.0 7 HAVG 

cee 0.5D0XBETAXEQEXSINCXM(C3) )XSINCXMC3) ) XCRE+HP ) 
=CxC 

22 = BETA ® XMC1) &® XMC2) 

ALP = 0.00 

MO = 1 

NN = 8 z 

CALL MMBSIR(ZZ,ALP,NN,MO, BF,NZ) 

IQ = BFC(1) 

lig= BEC) 

TZ2=—Des 32 

V5 = 8r¢9) 

IG = BFC5) 

I5 = BFC6) 

I6 = BFC7) 

ly = SB ELS) 

FA = FRA 

FBR = FXB 

E = XM(2) 

C2W = COSC2.*XM(4)) 

S2W = SINC2. XXMC4)) 

CG4W = COSCG.*XMC4)) 

SGW = SINCG.*XMCG)) 


tees =-DELXCSQRT CBOXXMC 1) )) XRHOOXEXPC-ZZ-CXC2W) *C C1. +0. 25D0% 

CXC )X 

1CIO+EXI1 )+FAXCC1.+0.25D0XCC)¥CI1+EXI2) +0. 5DOXCKC (2. XI2-EX(I1-3.% 
213) D+FARCI1+13+2.XEXIG) )XC2W-0. 5SDOXCXFRXCCI1-I13)+2.XEXCI 2-19) )% 

Sy ScWet be ZO DOXCKCK(CZ A1GaE eC 5 415-5. XIS)PFARCCISTISIFEXCS.XI6-12 
G4 JIIRCGN + O.L25DOXCKCKFBXCCIS-IS)+EXCI 2-4. X14+3 .¥16)) XS4W) 


XT1 = DEL XCSQRTCBOXXMC1)))XRHOO 
XT2 = DELMG 7 XTl 

Ato. = EXPC-2Z)4C10 + ExT1) 

Kus = XN2e7 ATS 


DELMG = 1090. * DELMG 7 9.8D0 

PRINT ®, ZZ, RHOO, XT1l, XT2, XT3, XT4, DELMG , DEL 
RETURN 

END 


HK HK KK MKA NRK KKK KKK RKHK KX KKK KK KKK 
PESTS E SES SESS SSS SESS SE SST SSE ES F | 


*X %X 
XX SUBROUTINE JAC60 xX 
x | xX 


HRN RKKKKRKRKKKRKRKKRKRKKRKRRKKRARKKKR 
RRR KKKRKRKRKRRRRKKRKKKKRKKRKRKRKRRR 


SUBROUTINE JAC60 (D,DS,A,AS,BLA,Z,F107,RHO,SH ) 

Bee= Seaxs19159D0 7 180. 

XP = -16.021D0 - 0.001985D0%Z + 6.363D0XEXPC -0.0026D0Z) 
RO = 10.%XXP 

CSPSI = SINCD)XSINCDS) + COSCD)XCOSCDS)XCOSCA-AS-BL) 

Ree (OpoweoxCl., + CSPSI))XX3. 

Rice fo fen tODOXCEXPCO.0055D0xZ) - 1.900) XF 
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mK KOM OK MK OK OK OK OK OK OK 


RHO = 0.01D0 * F107 * RO * RHOL X* (10.12) 


DRODZ = 2.3026D0*XROXC-.001985D0 - .016544D0XEXP(-.0026D0*Z) ) 
DR1IDZ = .OOLOG5DOXFXEXPC .0055D0%Z) 

SH = -CROXRHO1)/ CROXDRIDZ+RHOIXDRODZ) 

tl 

ND 


KK RKKAKKARKRAAKAKAKK KKK KKK AK AKKAKKKKKR 
HH HH HHH HHH HHH HH HH HH HHH HK HH HK 


XX %X 
KK SUBROUTINE THRST xX 
xX xX 


RH HK HHH ARK KHH KKK KKK KK KK KARA KK RKKK 
LoS SS SSeS 52 SS F225 2 $5535 555 25%: 


SUBROUTINE THRSTC XM, DXM, IRC, K, PA, DWOA ) 

DIMENSION XMC6) , DXMC6), IRVC1O), DAC10) 

COMMON 7 OADJ 7 NOA, IRV, DA, SPI, WT, PSENS 

DATA BO 7 398603.254D07, P1I/73.14159D0/ 

DO LO 1 = 7156 

DXMCI) = 0.0 

CONTINUE 

DWOA = 0.0 

IF € IRC .EQ. 0) RETURN 

THT ee + 2. ¥XMC2)¥*SINCXMC6)) + €5.74.)XXMC2) XXMC2)XSINC2Z. XXM 


R = XMC1)¥C1. -— XMCZ)RXMC2])IZC1. + XMC2)*COSCTHT)) 
V2 = BO X* €C€2.7R) = C1./7XMC1))) 

V = SQRTCV2) 

DVI = DACK)¥BO/€2.¥XMC1) ¥XMC1) ¥V) 

DXMC1)= DACK) 

DXMC2) = 2.*¥CCCOSCTHT)+XMC2) )/V) XDVI 

DXMCG) = 2.*¥SINCTHT)XDVI/CXMC 2) ¥V) 

DXM(6) = €2.*PI/PAJ“-C1l./CXMC1)¥CSQRTC1.-XMC2) XXMC 2) )) XV) D2. 
1 SINCTHT)¥C CXMC1) ¥O1. -XM02) XXMO2))/7XMC2)) + R¥XMC2) DVI 

PA = 3.*¥PAXXMC1)¥CV/B0)XDVI 

DNOA = CWNT/CSPI*9 .8D0))¥1000.xDVI 

IF CPSENS .EQ. 0) THEN 

PRINT ee SD Cee 

i DNOA 

END IF 


RETURI 
END 
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APPENDIX B 


SAMPLE INPUT AND OUTPUT 





05000,001,0,1,0,025 
2-2, -G000500, 1.0; 250%. oo ,09100. 
POOe, WwOOs 2.0 ,935222.758 
1 1094/7 U 90070 78 060 a 0 POO rG 791.1015 0 O ++ | 
Celg7s? 446275757 160//5 150.0000°9095.0000 200.0000 .O100000 095.0000 ++ | 
See toeeeoso7o4 6000S S92 F.@Oll2 =-3.81308 -0.23959E-00 0 ++ | 
4 10947 000000000-0 0 0 0 0 ++ | 
a 105947 Ole 0oo00000 -—701.8652 E-00 sna 
CCCO X.X XXXXXK . KKXXXXXXX XXXXXX .XXXKXXXXX 
CD 1.0 000001.000000000 000003.500000000 
A 0.0 000000.000010000 000000.001000000 
D 0.0 000000.100000000 000010.000000000 
SPI 0.0 000200.000000000 003000 .000000000 
W 0.0 000005.000000000 000200.000000000 
WT 0.0 000100.000000000 100000.000000000 
F107 0.0 000050.000000000 000300.000000000 
FBAR 0.0 000050.000000000 000300.000000000 
AKP 0.0 000000.000000000 000007.000000000 
TVE 0.0 038422.000000000 039422.000000000 
UJD 0.0 043223.000000000 044223.000000000 
cETU 0.0 000000.009000000 000179.000000000 
HO 0.0 000000.00000C000 000360.000000000 
GO 0.0 000000.000000000 060360.000000000 
ES 0.0 000000.010000000 000000.100000000 
me 0.0 009022.500000000 000157 .500000000 
RND 0.0 000000.000010000 000000.000000000 
ESD 0.0 000009.000000900 000000.000000000 
AO 0.0 000001 .0259490000 000001.168470000 
ADOT 0.0 000000.000090000 000000.000000000 
1000, 0.1 
2000, 0.3 
3000, 0.5 


KK KRAKA KKK RAK KK KKK KKK RRR RK KK RK RRR KKK RK KKK KK KKK KK KK KK KKK KK KK KKK RK AR ARK RKKRKKR 


x 5 3 
x OMEY THE ABOVE IS ACTUAL ENPUT FILE DATA. AN INPUT ELEMENT x 
AS OESCRIPTIVM BS PRESENTED WEXT TO AID THE USER IN PLACING VALUES ¥ 
x x 


Pe PS SSSSSESSSS SESS SSSESS SESE SS SESTESE SESE SSESESESTESSESSESE SESS SESESS SSE SESE? SF 3 


INPUT FILE FOR PROPELLANT LONGEVITY MODEL CVARIABLE ORGANIZATION) 


Pepne> IETRZ, NOA, PSENS, IET, PSIZE 

CD, Ay D, Sei W, WT, 7 

FlO7, FREAK, AKP, TVE, 
1 SAIP LINE ONE NAYSPASUR FIVE CARD DATA 
Coote, xe, 1X, RO, 1X, GO, 1X, ES, 1X, XI 

a eu, AND, ton ESD 
4 SKIP LINE FOUR NAVSPASUR FIVE CARD DATA 

5 OX; aluiy UX, ADOT 

TRVC1), DACI) 

meted, DAG?) NOTE: WARNING OMIT ALL ORBIT ADJUST DATA IF 


! NO ORBIT ADJUSTS ARE SCHEDULED (NOA=0) 
| 

IRV(X), DACK) | 

SENSITIVITY ANALYSIS TABLE 


CHAR SNSAR SNSMIN SNSMAX 
AS meee FEGaiG Elen 6 
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=_— ep op @p —* oe oF om 6 22 a 62 ee a oP ee a ee ow op Se 2 a 6 ow ee om 2 6 Se aS ew ow ow © ow Se © Oe SS Sw Se Se es Ss Ss 


—_— es =e oe oe oe a oe a 2D 2 6 22 22 6 62 6 2 oe 2 oe oe ew oe aS 2 a SS > 6D GD 6D GD OD CD SOD CD OD SD SD OD OD Se GD OD a oF 


ATMOSPHERIC INPUT PARAMETERS 


D CRELATIVE ATMOSPHERIC ROTATION): = 1.0000 

F1O7 CDECIMETRIC SOLAR FLUX) = A aban: 

FBAR (90 DAY AVG. F107) = aie 

AKP CGECMAGNETIC INDEX) = 2. 

TVE CTIME OF VERNAL EQUINOX PASSAGE) = 43222.738200000 MODIFIED JULIAN DAYS 


SATELLITE PHYSICAL PARAMETERS 


CD CDRAG COEFFICIENT) = 2.20 UNITLESS 
A (CCROSSECTIONAL AREA) = 0.000050000 SQUARE KILOMETERS 


SPI CMOTOR SPECIFIC TMPUESED R= 230.00 aa 
W CINITIAL FUEL MASS) = 10.00 KILOGRAM 
WT CINITIAL SATELLITE MASS) = 9100.00 XILOGRAMS 


NAVSPASUR ORBITAL ELEMENTS 


UJD CEPOCH) = 44619.98716775 MEAN JULIAN DAYS 

ETU (MEAN ANOMALY) = 150.0000 DEGREES 

HO (MEAN RIGHT ASCENSION OF THE ASCENDING NODE) = 95.0000 DEGREES 

GO (MEAN ARGUMENT OF PERIGEE) = 200.0000 DEGREES 

XI (MEAN INCLINATION) = 95.0000 DEGREES 

RND CRATE OF CHANGE OF MEAN MOTION) = 0.000783912 REVOLUTTIONS PER DAY SQUARED 


ES (MEAN ECCENTRICITY) = 0.0100 UNITLESS 

ESD CECCENTRICITY DECAY RATE) = -0.239590000 I/SECONDS 
AO CKAULA SEMI-MAJOR AXIS) = 1.06000 MEAN EARTH RADII 
ADOT CSEMI-MAJOR AXIS DECAY RATE) = -701.865200000 


CONTROL INPUT PARAMETERS 


MAXIMUM REVOLUTIONS PER ITERATION 5000 
NUMBER OF SCKEDULED ORBIT ADJUSTS TO SEMI-MAJOR AXIS 0 


SENSITIVITY ANALYSIS FOR 
LOW EARTH ORBIT SATELLITES 
PROPELLENT LONGEVITY 
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TABLE OF INPUT PARAMETERS SCANNED FROM MINIMUM TO 
MAXIMUM WITH THE INDICATED INCREMENT PER RUN. 

THE GRANULARITY IS 1/7 25. 

ALL OTHER INPUT PARAMETERS ARE RETURNED TO THEIR 
INITIAL VALUES AT EACH INCREMENT. 


INPUT PARAMETER MINIMUM MAXIMUM INCREMENT 
CD 1.000000000 3.500000000 0.100000000 


SENSITIVITY SCAN RESULTS FOR ABOVE SPECIFICS 


REV NUMBER TIMECDAYS) FRACTION OF SCAN 
495 Slo6 _ 1.000000000 
448 237 1.100000000 
409 26.08 1.200000000 
376 eon 77 1,.300000000 
348 Zone r.400000000 
324 20.65 1.500000000 
303 19230 1.600000000 
285 F3.15 1.700000000 
268 be ge 1s 1.800000000 
254 Mee 17 1.900000000 
241 15.354 2.000000000 
cco 14.57 2.100000000 
2h6 £35.87 2.200000000 
208 PS25 2.300000000 
200 12512 2.400000000 
Meal 12.14 2.500000000 
184 Tio 2.600000000 
ar 1 ras, 2.700000000 
179 10.30 2.800009000 
1S QRZ 2.900000000 
59 1 Oen0 3.0900000000 
154 2.76 3.100000000 
149 9.56 3.200000000 
144 9.14 3.300000000 
140 8.88 3.400000000 
136 8.63 3.500000000 


Pony 


05000,001,0,1,0,025 
2.2, .0000500, 1.0, 250.5001. 09.00- 
100., 100., 2.0 ,S322257o62 


1 10947 VU 90070 78 060 A 0 10 01S7 91-1015 0 0 ++ | 
2 10947 44619.98716775 150.0000 095.0000 200.0000 .0100000 095.0000 ++ | 
5 10967 15.71356734 .000763912 1.00112 °-35 °s15065-0 23727 e-00 QO ++ | 
4 10947 000000000-0 0 0 0 0 ++ | 
5 10947 01.06000000 -701.8652 E-00 ++ | 
CCCC X.X XXXXXX. XXXXXXXXX XXXXXK . XXXXXXXXX 
CD 0.0 000001.000000000 000003.500000000 
A 0.0 000000.000010000 000000.001000000 
D 0.0 000000.100000000 000010.000000000 
SPI 0.0 9000200.000000000 003000.000000000 
W 0.0 000005.000000000 000200 .000000000 
WT 0.0 000100.000000000 100000.000000000 
F107 0.0 000050.000000000 000300.000000000 
FBAR 0.0 990050.000000000 000300 .000000000 
AKP 0.0 000000.000000000 000007 .000000000 
TVE 0.0 038422.000000000 0394922.000000000 
UJD 0.0 043223.000000000 044223.000000000 
ETU 0.0 000000 .000000000 000179.000000000 
HO 0.0 000000.000000000 000360 .000000000 
GO 0.0 000000.000000000 . 000360.000000000 
ES 0.0 000000.010000000 000000.100000000 
XI 0.0 000022.500000000 000157 .500000000 
RND 0.0 000000.000010000 000000 .Q000000000 
ESD 0.0 000000.000000000 000000.000000000 
AQ 1.0 000001.025940000 000001.168470000 
ADOT 0.0 000000.000000000 000000.000000000 
1000, Q.1 
2000, 0.3 
3000, 0.5 


HHH HR HH HH RH KKH HK HH HHH HHH HK HHH HHH HHH HH HHH HR K RH HHH HR HH KK KKK KKK KKK KKK KKK KR KKKRK 


% % 
* ONLY THE ABOVE IS ACTUAL INPUT FILE DATA. AN INPUT ELEMENT x 
* DESCRIPTION IS PRESENTED NEXT TO AID THE USER IN PLACING VALUES x 
¥% % 


HEHE DE HE DE HE HE HE HE HE DE HE HE HEHE HE I HE IE HE IE HE HE IE 3 HE HE HE HE HED IE HE DE IE DE HE HE IE HE HE HE I HE HE HEE HE DK HE HE DE DE HE DE DE HE DE IE HE HE HE DE HE HE IE DE DE HE HE EK 
INPUT FILE FOR PROPELLANT LONGEVITY MODEL CVARIABLE ORGANIZATION) 


TETRIG SLEt RZ NORA, ESENS LE, PSIZe 

CD, A, D, SPI, W, WT, 

F107, FBAR, AKP, TVE, 

1 SKIP LINE ONE NAVSPASUR FIVE CARD DATA 
2 8X, VUJD, IX, EVU, 1%, HO, 1X, GO, Pees. 1x, 
- 3 20X77 eel Bie sx ESD 
G SKIP LINE FOUR NAVSPASUR FIVE CARD DATA 
5 8X, AO, 1X, ADOT 


IRVO1), DAL) 
IRVO2), DAC2) NOTE: WARNING OMIT ALL ORBIT ADJUST DATA IF 
33, DAC33 NO ORBIT ADJUSTS ARE SCHEDULED (CNOA=0} 


| 
| 
IRVOK), DACK) | 
SENSITIVITY ANALYSIS TABLE 


CHAR SNSAR SNSMIN SNSMAX 
AG gol! F16.8 F16.8 


FE DE HE FE IE IE IESE HEE EE DE HE IE DE HE FE HE IE EE IE DE IE IE IE IE IE IE IE IE HE HE HE HE HE HE HE EE HE DE HE IE HE HEE HE EE HEE EE EK HR RR RK NRK 
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aS ee SS SS SP SP SP SP SP SP SB PSP VOWS PWMPS#8 SP SBP SP SF SF S28 SF SP SF SF FOS OP SS 2 @ @S @S wo @=S = = 


ATMOSPHERIC INPUT PARAMETERS 


D CRELATIVE ATMOSPHERIC ROTATION) = 1.0000 

F197 CDECIMETRIC SOLAR FLUX) = 100.00 

FeAR (90 DAY AVG. F107) = 100.00 

AKP CGEOMAGNETIC INDEX) = 2.000 

TVE (TIME OF VERNAL EQUINOX PASSAGE) = 43222.738200000 MODIFIED JULIAN DAYS 


SATELLITE PHYSICAL PARAMETERS 


CD (DRAG COEFFICIENT) = 2.20 UNITLESS 
A CCROSSECTIONAL AREA) = 0.000050000 SQUARE KILOMETERS 


Sal GrOUGcworecirlC IMPULSE) = 220, 0UR SECONDS 
WeCINITIAL FUEL MASS) = 1.00 KILOGRAMS 
WT CINITIAL SATELLITE MASS) = 9100.00 KILOGRAMS 


NAVSPASUR ORBITAL ELEMENTS 


UJD CEPOCH) = 44619.98716775 MEAN JULIAN DAYS 

ETU (MEAN ANOMALY) = 150.0000 DEGREES 

HO (MEAN RIGHT ASCENSION OF THE ASCENDING NODE) = 95.0000 DEGREES 

GO (MEAN ARGUMENT OF PERIGEE) = 200.0000 DEGREES 

XI (MEAN INCLINATION) = 95.0000 DEGREES 

RND CRATE OF CHANGE OF MEAN MOTION) = 0.000783912 REVOLUTTIONS PER DAY SQUARED 


Ear viEANeceeENTRICITY) = 0.0100 UNITLESS 
ESDeCECCCNMRICITY DEGAY RATE) = -0.239590000 1/SECONDS 
AO (KAULA SEMI-MAJOR AXIS) = 1.06000 MEAN EARTH RADII 
ADOT (CSEMI-MAJOR AXIS DECAY RATE) = -701 .3865200000 


CONTROL INPUT PARAMETERS 


MAXIMUM REVOLUTIONS PER ITERATION 5000 
MUrseER OF SCHEDULED ORBIT ADJUSTS TO SEMI-MAJOR AXIS 0 


meee eee O22 SO SP 2 SPP SPV SP SP SP PSP SPP PSP SP 2 SP SP SP SP PSPs ae Sow Pes V2 V2 2 SS Ss = == | 


SENSITIVITY ANALYSIS FOR 
LOW EARTH ORBIT SATELLITES 
PROPELLENT LONGEVITY 
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-_ en oe oP ee Se ee 2 es? 6 Se 2 ee ee ee ee Se Se Se Se ee ee a ee eS a 6 ae oo Se So 6 Se Se Se ow ee Ge ee oe ee ee ee oe a ee = == 


TABLE OF INPUT PARAMETERS SCANNED FROM MINIMUM TO 
MAXIMUM WITH THE INDICATED INCREMENT PER RUN. 
THE GRANULARITY Y ise 7eo. 
ALL OTHER INPUT PARAMETERS ARE RETURNED TO THEIR 
INITIAL VALUES AT EACH INCREMENT. 
INPUT PARAMETER MINIMUM MAXIMUM INCREMENT 


AO 1.025940000 1.168470000 0.005701200 


@2[2@ 2S 2 SP SP SP SP PSPSPS OPPs SP SPP SSeS OPPS SPP SFP SP SF SFP SP 2 S&S SS SS SS & = = = 


SENSITIVITY SCAN RESULTS FOR ABOVE SPECIFICS 


REV NUMBER TIMECDAYS) FRACTION OF SCAN 
l 0.00 1.025940000 

i 0.00 1.031641200 

1 0.00 1.037 342400 

> O12 1.043043600 

6 de Sih 1.0438744800 

YA 0.70 1.054446000 
22 1. 34 1.060147200 
40 Pape | 1.065848400 
69 4.42 1.071549600 
TTS 7.497 1.077250800 
186 Erte | 1 0629452 0:00 
298 19.76 1.088653200 
473 Sl 365 1.094354400 
747 50.42 1.100055600 
1149 ere, 1.105756300 
1652 Tis. 55 1.1114538000 
2055 140.55 per Oh Lays) raga) 1 
2270 138.15 1.122860400 
2514 Ss he 1.1285616C0 
2942 208-15 1.154262300 
3349 250-72 1.1399640900 
3645 262.09 1.145665200 
4026 Zo oe 1 .152266900 
4373 Sl2.60 1.157067600 
4555 334.53 1.162768300 
4712 348.62 1.163470000 
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05000,001,0,1,0,025 
Even OUNe S00, 1.0, 250. PrO? EGO. 
E00. ees, 2-0. ,495222.73382 
i 10947 U 90070 78 060 Poors; 9!) 101s 0 0 | 
2 10947 44619.938716775 150.0000 Bren: 200.0000 .0100000 095.0000 | 
3 10947 15.71356734 ponies SU Zene a5 bo0S =0.25959E-00 | 
4 10947 000000000-0 0 | 
5 10947 01.06000000 -701.8652 e- 00 | 
CCCOC X.X XXXXXXK.XXXKXXKXXX XXXXXX .XXXXXXXXX 
CD 0.9 000001.000000000 000003.500000000 
A 0.0 000000.000010000 000000.001000000 
D 0.0 000000.100000000 000010.000000000 
SPI 90.0 000200.000000000 003000.000000000 
W 0.0 000005.000000000 000200.000000000 
WT 0.0 090100.000000000 100000.000000000 
F107 1.9 000050.000000000 000300.000000000 
FBAR 0.0 000050.000000000 000300.000000000 
AKP 0.0 000090.060000000 000007.000000000 
TYE 0.0 038422.0000000C60 032422.000000000 
UJD 0.0 045223.000000000 044¢223.000000000 
ETU 0.0 000000.000000000 000179.000000000 
HO 0.9 000000.000000000 000360.000000000 
GO 0.9 000000.000000090 000360.000000000 
ES 0.0 000000.010000000 000000.100000000 
XI 0.0 000022.500000000 000157 .500000000 
RND 90.0 000000.000010000 000000 .000000000 
ESD 0.0 000000.000000000 000000.000000000 
AO 0.0 000001.025940000 000001.168470000 
ADOT 0.0 000000.000000000 000000 .000000000 
1000, 0.1 
2000, 0.3 
3000, 0.5 


HH HH HHH HEH HEH HEHE HE HEHEHE HE HK HE HH HE HE HE HE HE HEH HEE HEH HK HK HK HK HH HH HK HH KK KK RR KRAARA RRR RARRRARRARKKRARKRK - 


x x 
* ONLY THE ABOVE IS ACTUAL INPUT FILE DATA. AN INPUT ELEMENT 
* DESCRIPTION IS PRESENTED NEXT TO AID THE USER IN PLACING VALUES X 
X x 


HEH HEHEHE HEH HEHE HE HE HEHE HE HE HE HE HE HEH HE HE HE HE HEH HE HEHE HE HE HE HE HE HE IE HE HE HE HE HEH HE HE HE HE HE HE HE HEHE HE DE HE HE DE HE HE DE HE DE DE HE HE HE EK EK KH 
INPUT FILE FOR PROPELLANT LONGEVITY MODEL (VARIABLE ORGANIZATION) 
DehRis LEGR<, NOR, PSENS, LET, PSIZE 


CD: A, D, SPI, W, ? 
F107, FBAR, AKP, TVE, 


1 SKIP LINE ONE NAVSPASUR FIVE CARD DATA ’ 
2 3X, UJD, 1X, ETU, 1X, HO, 1X, GO, 1X, ES, 1X, XI 

3 20X, RND, 19X, ESD 

G@ SKIP LINE FOUR NAVSPASUR FIVE CARD DATA 

5 8X, AO, 1X, ADOT 

IRVC1), DACI) 

IRVO2Z), DAC2) NOTE: WARNING OMIT ALL ORBIT ADJUST DATA IF 

PRvCS), NO ORBIT ADJUSTS ARE SCHEDULED CNOA=0) 


| 
| 
| 
| 
IRVOK), DACK) | 
SENSITIVITY ANALYSIS TABLE 


CHAR SNSAR SNSMIN SNSMAX 
AG me) a Fl6.8 F16.8 


HH HK HHH ARR HK RR KK RR RHR RR K KR RRR HHH HK HHH HHH HHH HK HEHE HEH EH HK KR KK HK RK KH KKK KR ARARRRKRKR 
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man as am am am am ae am am ae a 22 22 GF a 22 a> 22 22 22 22 2E OE ae oe 2D oe ae GE GE GE GE GP GP GP GP GE GP OP 4D om 6S 2 SF ae ae = 


ATMOSPHERIC INPUT PARAMETERS 


D CRELATIVE ATMOSPHERIC ROTATION) = 1.0000 

F197 CDECIMETRIC SOLAR FLUX) = 100.00 

FBAR €90 DAY AVG. F197) = 100.00 

AKP CGEOMAGNETIC INDEX) = 2.000 

TVE CTIME OF VERNAL EQUINOX PASSAGE) = 43222.738200000 MODIFIED JULIAN DAYS 


SATELLITE PHYSICAL PARAMETERS 


CD (DRAG COEFFICIENT) = 2.50 UNITLESS 
A CCROSSECTIONAL AREA) = 0.000050000 SQUARE KILOMETERS 


SPI CMOTOR SPECITETC IMPULSE =. 230.00 SECONDS . 
W CINITIAL FUEL MASS) = 10.00 KILOGRAMS 
WT CINITIAL SATELLITE MASS) = 9100.00 KILOGRAMS 


NAVSPASUR ORBITAL ELEMENTS 


UJD CEPOCH) = 44619.98716775 MEAN JULIAN DAYS 

ETU (MEAN ANOMALY) = 150.0000 DEGREES ; 

HO (MEAN RIGHT ASCENSION OF THE ASCENDING NODE) = 95.0000 DEGREES 

GO (MEAN ARGUMENT OF PERIGEE) = 200.0000 DEGREES 

XI (MEAN INCLINATION) = 95.0000 DEGREES 

RND CRATE OF CHANGE OF MEAN MOTION) = 0.000783912 REVOLUTTIONS PER DAY SQUARED 


ES (MEAN ECCENTRICITY) = 0.0100 UN IahiskSsS 

ESD CECCENTRICITY DECAY RATE) = -0.239590000 1/SECONDS 
AQ (KAULA SEMI-MAJOR AXIS) = 1.06000 MEAN EARTH RADII 
ADOT CSEMI-MAJOR AXIS DECAY RATE) = -701.865200000 


CONTROL INPUT PARAMETERS y 


MAXIMUM REVOLUTIONS PER ITERATION 5000 
NUMBER OF SCHEDULED ORBIT ADJUSTS TO SEMI-MAJOR AXIS 0 


2 a am am am 2 GP GP GP OE OP OE SE OE SESE SST SI SP PBS VSS SF SF SSF SBSBSB VS SBSOQSSB SB SPSS SS SS |] = 


SENSITIVITY ANALYSIS FOR 
LON EARTH ORBIT SATELLITES 
PROPELLENT LONGEVITY 
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TABLE OF INPUT PARAMETERS SCANNED FROM MINIMUM TO 
MAXIMUM WITH THE INDICATED INCREMENT PER RUN. 
THE GRANULARITY IS 1/7 25. 
ALL OTHER INPUT PARAMETERS ARE RETURNED TO THEIR 
INITIAL VALUES AT EACH INCREMENT. 
INPUT PARAMETER MINIMUM MAXIMUM INCREMENT 


ars 50.000000000 300.000000000 10.000000000 


SENSITIVITY SCAN RESULTS FOR ABOVE SPECIFICS 


REV NUMBER TIMECDAYS) FRACTION OF SCAN 
392 . 29.499 50.000000000 
324 20-00 60.000000000 
276 17.58 70.000000000 
241 15. 54 80.000000000 
Fags’ Lo. lS 90.000000000 
191 T2. 14 100.000000000 
174 11.06 110.000000000 
159 ro-10 120.000000000 
146 ed 130.000000000 
JESS 8.63 140.000000000 
127 S.0> 150.000000000 
119 7.594 160.000000000 
1ll2 7 170.000000000 
105 6.65 180.000000000 
100 G5 190.000000000 

25 @. 01 200.000000000 
50 5.69 210.000000000 
86 5.43 220.000000000 
82 els 230.00000000CG 
79 4.99 240.000000000 
73 4.79 250.000000000 
es 4.60 260.000000000 
70 4.41 270.000000000 
67 4.22 280.000000000 
65 4.09 290.000000000 
63 96 300.000000000 
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ATMOSPHERIC INPUT PARAMETERS 


‘D CRELATIVE ATMOSPHERIC ROTATION) = 1.0000 

F107 CDECIMETRIC SOLAR FLUX) = 100.00 

FBAR (90 DAY AVG. F107) = 100.00 

AKP CGEQMAGNETIC INDEX) = 2.000 

TVE CTIME OF VERNAL EQUINOX PASSAGE) = 43222.738200000 MODIFIED JULIAN DAYS 


SATELLITE PHYSICAL PARAMETERS 


CD (DRAG COEFFICIENT) = 2.20 UNITLESS 
A CCROSSECTIONAL AREA)d = 0.000050000 SQUARE KILOMETERS 


SPI CMOTOR  SPEGTETC Ciirurse) = 230.00 SECONDS 
W CINITIAL FUEL MASS) = 150.00 KILOGRAMS 
WT CINITIAL SATELLITE MASS) = 9100.00 KILOGRAMS 


NAVSPASUR ORBITAL ELEMENTS 


UJD CEPOCH) = 44619.98716775 MEAN JULIAN DAYS 

ETU CMEAN ANOMALY) = 150.0000 DEGREES 

HO (MEAN RIGHT ASCENSION OF THE ASCENDING NODE) = 95.0000 DEGREES 

GO (MEAN ARGUMENT OF PERIGEE) = 200.0000 DEGREES 

XI (MEAN INCLINATION) = 95.0000 DEGREES 

RND CRATE OF CHANGE OF MEAN MOTION) = 0.090783912 REVOLUTTIONS PER DAY SQUARED 


ES (MEAN ECCENTRICITY) = 0.0100 UNITLESS 

ESD CECCENTRICITY DECAY RATE) = -0.239590000 I“SECONDS 
AQ CKAULA SEMI-MAJOR AXIS) = 1.06000 MEAN EARTH RADII 
ADOT CSEMI-MAJOR AXIS DECAY RATE) = -701.865200000 


CONTROL INPUT PARAMETERS 


MAXIMUM REVOLUTIONS PER ITERATION 5000 
NUMBER OF SCHEDULED ORBIT ADJUSTS TO SEMI-MAJOR AXIS 0 


PROPELLANT LIFE ESTIMATOR 
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GDA= 0.1) 00G8E-03eSQ KM 
SeeClEle thr vesSeso250.00 SEC 
INIVIAL WEIGHT OF PROPELLANT= 150.00 KG 


F107= 100.00 
FBAR= 100.00 
TIME OF VERNAL EQUINOX= 0.43222738200000E+05 


MAXIMUM NUMBER OF ORBITAL REVOLUTIONS= 5000 


DATA PRINTED EVERY 100 ORBITAL REVOLUTIONS 


KOZAIT MEAN ELEMENT SET 


SEMIMAJOR AXIS= 0.67562054107790E+04 KM 
EGRENTRICTIY= 0.96561152904615E-02 
INCLINATION= 0.94999964849617E+02 DEG 
ARGUMENT OF PERIGEE= 0.19333742463643E+03 DEG 
ASCENDING NODE= 0.95000149668470E+02 DEG 

MEAN ANOMALY= 9.15661152641309E+03 DEG 


REV NUMBER TIMECMJD) TIMECDAYS) PROPELLANT REMAININGCKG) 
100 44626 . 32 6.55 145.28 
200 449632.71 Vea 2 140.69 
300 44639.10 Model 136.021 
400 44645 .49 25.50 LSS 
500 44651.88 Sl 0 127.54 
600 44658.27 38.29 Jig 
700 44664.67 44.68 LEGS 
800 44671.06 So? Oey) 
900 44677 .45 57.46 NE Ore 72 

10090 44683.84 Gai.o9 108.907 
1100 44690.23 7 Oe25 104.42 
1200 44696 .63 76.64 TUDE 7 2 
1300 44703 .02 69.03 96.94 
1400 44709.41 89.42 SNES Us, 
1500 44715.80 Sra 2 Sa- 96 
1600 44722.19 hOZ. 21 84.60 
1700 44728.59 108.60 he ee, A 
1360 44734.98 a9 9 74.98 
1960 44741.37 hess 69725 
2000 44747 .76 ier 69.97% 
2100 44754.15 TsS.17 GOn 27 
2290 44760.55 149.56 ID.92 
2200 44766 .94 146.95 BRE 
2400 44773.33 153.34 48.18 
2500 44779 .72 LoS 44.59 
2600 44786 .11 66.15 41.07 
2700 44792.51 Tae oe of oS 
2800 44798 .90 18.91 34.10 
2900 44805.29 ToS 230 50.53 
3000 44811.68 191.69 Zora! S 
3100 44818 .07 198.09 22.74 
3200 44324 .47 204.43 oor 
3300 44830 .86 210.87 SISK, 
3400 44357 .25 2726 8.68 
3500 44843 .64 EoougGS ST gene 
3600 44850.03 Ca0'.05 Go 
3600 230.65 0.000000000 
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